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Abstract 


In  the  framework  of  the  proposed  "continuous  approach”  to  constrai¬ 
ned  optimization  problems,  we  describe  two  new  solution  methods  which  re¬ 
sulted  from  the  research.  The  first  is  a  continuous  "inexact"  method  for  sol¬ 
ving  systems  of  nonlinear  equations  and  complementarity  problems  (along 
the  lines  of  the  DAFNE  Method),  and  the  second  is  a  continuous  method 
for  solving  the  linear  programming  problems  (along  the  lines  of  Karmarkar’s 
method)  which  is  shown  to  be  quadratically  convergent. 

Some  numerical  experience  on  a  number  of  test  problems  is  reported. 
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1.  Introduction 

This  is  the  final  report  on  the  work  performed  from  September  1986 
to  December  1991,  under  contract  n.  DAJA  45-86-C-0028  awarded  to  the 
University  of  Rome  "La  Sapienza"  on  the  research  project  "Numerical  Opti¬ 
mization",  by  the  principal  investigator  Francesco  Zirilli  and  his  co-workers. 

The  objective  of  the  research  is  described  in  par.  2,  the  results  of  the 
research  are  described  in  par.  3,  and  some  conclusions  are  in  par.  4. 


2.  Objective  of  the  research 

The  subject  of  the  research  was  the  field  of  those  problems  in  con¬ 
strained  optimization  which,  starting  from  the  linear  programming  problem, 
can  be  formulated,  with  growing  degree  of  generalization,  first  as  linear 
complementarity  problems  and  second  as  nonlinear  complementarity  pro¬ 
blems. 

The  objective  of  the  research  was  to  attack  the  above  problems  by 
means  of  the  so-called  "continuous  approach"  to  optimization  (as  opposed  to 
the  so-called  "pivotal"  methods,  such  as  the  simplex  method  for  linear  pro¬ 
gramming),  with  special  consideration  for  the  interesting  cases  of  non-con- 
vex  or  ill-conditioned  problems,  and  problems  with  a  very  large  number  of 
variables;  and  in  particular  the  objective  was  to  investigate  the  possibility  of 
applying  to  the  above  problems,  suitably  transformed  into  no  ilinear  equa¬ 
tion  problems,  the  methods  developed  by  the  principal  investigator  and  his 
co-workers  for  solving  nonlinear  equations  and  global  optimization  pro¬ 
blems,  based  on  the  numerical  integration  of  suitable  ordinary  or  stochastic 
differential  equations  (refs.  [1]  to  [4]). 


3.  Results  of  the  research 


During  the  development  of  the  research  the  complementarity  pro- 
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blems  proved  to  be  much  more  difficult  than  it  had  been  anticipated  and  it 
became  clear  that  the  original  plans  where  somehow  too  ambitious. 

The  final  outcome  of  the  research,  if  judged  against  the  original 
plans,  is  therefore  admittedly  less  satisfactory  than  it  was  originally  hoped; 
nevertheless  a  number  of  interesting  results  have  been  obtained,  so  that  we 
feel  that  the  research  is  still  to  be  considered  at  least  partially  successful. 

The  main  results  of  the  research  are  contained  in  the  two  papers 
numbered  [5]  and  [6]  in  the  list  of  references,  which  are  described  in  the  fol¬ 
lowing  paragraphs  3.1  and  3.2,  and  enclosed  as  Appendix  1  and  Appendix  2. 

Report  on  some  work  performed  in  other  directions  along  the  lines  of 
the  original  research  plan,  together  with  some  related  results  of  auxiliary 
and  preliminary  nature,  were  described  in  the  Periodic  Technical  Reports; 
see  also  the  papers  numbered  (8]  and  [9]  in  the  list  of  references. 

The  research  has  also  stimulated  scientific  contacts  with  several  Ita¬ 
lian  and  foreign  scholars. 

The  above  results  have  been  disseminated  by  means  of  the  aforemen¬ 
tioned  papers  on  high-standard  academic  journals,  and  seminars  at  Accade- 
mia  dei  Lincei,  Rome  (ref.  [7],  which  originated  paper  [8]),  at  two  meetings 
of  CECAM,  Centre  Europeen  de  Calcul  Atomique  et  Mol^culaire,  the  first 
in  Ermelo  (The  Netherlands),  ref.  [10],  and  the  second  at  CECAM  main  of¬ 
fice  in  Orsay  (Paris),  France,  ref.  [11]. 


3.1.  The  first  paper 

The  first  paper  (ref.  [5],  and  Appendbc  1)  can  be  summarized  as  fol¬ 
lows. 

A  class  of  algorithms  is  developed  for  the  numerical  solution  of  non¬ 
linear  systems  of  equations  and  complementarity  problems,  based  on  the 
fact  that  the  solution  of  complementarity  problems  can  be  reduced  to  the  so¬ 
lution  of  systems  of  nonlinear  equations  by  means  of  a  transformation  first 
suggested  by  Mangasarian. 

The  method  is  "continuous"  since  it  looks  for  the  solution  of  the  non- 
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linear  system  by  following  the  numerical  solution  trajectories  of  a  suitable 
differential  equation,  as  in  previous  work  of  the  same  authors  such  as  the 
method  implemented  in  the  package  DAFNE,  described  in  Refs.  [1]  and  [2]. 

At  each  numerical  integration  step,  the  DAFNE  method  requires  the 
solution  of  an  NxN  system  of  linear  equations,  and  the  cost  of  solving  such  a 
system  when  a  large  numer  N  of  unknowns  is  involved  is  the  most  important 
part  of  the  computation. 

The  present  method  can  be  called  "inexact",  since  it  computes  only  an 
approximate  solution  of  the  above  linear  system,  by  means  of  a  conjugate- 
gradient  procedure  which  is  suitably  stopped  before  "convergence",  i.e,  after 
a  number  mC<N)  of  steps  depending  on  the  norm  of  the  residual.  For  these 
algorithms  local  convergence  and  Q-superlinear  rate  of  convergence  has 
been  proved.  The  algorithms  have  been  used  to  solve  three  complementarity 
problems  derived  from  variational  inequalities  of  mathematical  physics  very 
successfully.  The  complementarity  problems  considered  had  up  to  900  varia¬ 
bles. 


3.2.  The  second  paper 

The  second  paper  (ref.  [6]  and  Appendix  2)  can  be  summarized  as 
follows. 

The  paper  introduces  a  new  method  for  solving  the  linear  program¬ 
ming  problem,  i.e.  the  problem  of  minimizing  a  linear  cost  function  of  seve¬ 
ral  real  variables,  subject  to  linear  equality  and  inequality  constraints. 

Following  Karmarkar  [12],  the  paper  considers  the  problem  in  the 
"canonical"  form 

minimize  f(x)  =  jf^x 

subject  to 

Ax  =i) 

l^x  =  1 

Xj^O,  i  =  l, ....  n 
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where 

ic  = 

£  ~  (cjf—jCjj)  , 

are  real  column  vectors  with  n  elements,  A  is  a  real  m  x  n  matrix  of  rank  m, 
with  As=i),  n^2,m<  n,  and  without  loss  of  generality  the  objective  fun¬ 
ction  f(x)  may  be  "normalized",  i.e.  f(x*)  =  0  if  x*  is  a  solution  of  the  problem. 

In  this  paper  it  is  shown  that  Karmarkar’s  method  [12]  is  in  fact  equi¬ 
valent  to  applying,  to  a  suitable  initial  value  problem  for  a  system  of  ordina¬ 
ry  differential  equations,  the  numerical  integration  method  known  as  Euler’s 
method  with  variable  stepsize,  and  obtaining  the  problem  solution  x*  as  the 
limit,  as  t  goes  to  infinity,  of  the  numerically  computed  solution  x(t)  to  the 
initial  value  problem,  starting  from  the  initial  point  Xq  =(l/n)g. 

The  proposed  method  is  also  based  on  the  above  interpretation  of 
Karmarkar’s  method,  but  with  two  main  differences: 

1)  the  initial  value  problem  is  based  on  a  different  system  of  ordinary  diffe¬ 
rential  equations; 

2)  the  numerical  integration  method  is  a  linearly  implicit  A-stable  method 
with  variable  stepsize. 

The  resulting  algorithm  is  shown  to  be  quadratically  convergent. 

The  computational  cost  of  one  step  of  the  proposed  algorithm  is 
shown  to  be  of  the  same  order  of  one  step  of  Karmarkar’s  algorithm. 

While  one  step  of  the  classical  simplex  algorithm  [13]  for  linear  pro¬ 
gramming  is  much  cheaper,  it  may  be  expected  that  -  due  to  the  quadratic 
convergence  -  the  number  of  iterations  needed  to  solve  a  linear  program¬ 
ming  problem  to  a  given  accuracy,  should  be  approximately  independent  of 
the  problem  size  n. 

Some  numerical  results  are  also  reported,  which  appear  to  support 
such  expectation. 

The  algorithm  was  tested  on  ten  test  problems,  one  originating  from 
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the  operations  of  an  industrial  plant  in  central  Italy,  and  the  other  nine  pro¬ 
vided  by  the  System  Optimization  Laboratory  at  Stanford  University. 

The  results  are  reported  in  Table  1  where  n  is  the  number  of  varia¬ 
bles,  m  the  number  of  constraints,  k  is  the  index  of  the  first  step  that  verifies 
the  stopping  rule 


f(x)  j<  10"*  .  f(Xo) 


and  V  jj  is  the  corresponding  value  of  f(x). 

We  note  that  the  test  problems  with  n,m  _<  5  are  solved  in  about  ten 
steps,  and  that,  while  n  and  m  vary  by  an  order  of  magnitude,  the  number  k 
of  steps  needed  to  solve  the  problem  varies  only  by  a  factor  of  two. 

TABLE  1 


Test  problem 

m 

n 

k 

1.  ZIRl 

304 

543 

21 

2.41D-10 

2.  ADUTTLE 

57 

141 

21 

3.16D-09 

3.  AFIRO 

28 

54 

12 

1.52D-12 

4.  BEACONFD 

173 

298 

20 

3.91D-09 

5.  BLEND 

75 

117 

21 

1.47D-12 

6.  ISRAEL 

175 

319 

17 

1.94D-10 

7.  SC105 

106 

166 

13 

1.36D-11 

8.  SC50A 

51 

81 

14 

1.48D-14 

9.  SC50B 

51 

81 

11 

7.84D-10 

10.  SHARE2B 

97 

167 

21 

1.78D-10 

4.  Conclusions 


The  research  resulted  in  two  new  methods,  one  for  solving  comple¬ 
mentarity  problems,  and  the  other  for  solving  the  linear  programming  pro¬ 
blems,  both  based  on  the  so-called  "continuous"  approach  to  optimization. 
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Successful  solution  was  obtained  for  the  complementarity  problems 
on  test  problems  with  up  to  900  variables. 

However,  due  to  the  great  difficulty  of  complementarity  problems  -  in 
fact  much  greater  than  expected  -  the  hoped  for  attack  of  much  more  diffi¬ 
cult  problems  proved  to  be  unsuccessful. 

The  linear  programming  method  was  shown  to  be  quadratically  con¬ 
vergent,  and  was  successfully  tested  on  preliminary  test  problems  with  up  to 
about  300  variables  and  540  constraints. 
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APPENDIX  1 


Rendiconti  di  Matcmatica.,  Setie  VII 
Volume  9  ,  Roma  (1989),  521-543 


An  Inexact  Continuous  Method  for  the  Solution 
of  Large  Systems  of  Equations  and 
Complementarity  Problems 

F.  ALUFFI-PENTINI  -  V.  PARISI  -  F.  ZIRILLI^*) 


Dedkito  alh  memoria  di  Carlo  Cattaneo,  maestro  ed  amico 


RjaSSUNTO  -  Si  considero  un  nuovo  metodo  per  la  risoluzione  numerica  sia  di 
sistemi  di  equazioni  non  /mean'  5ia  di  problemi  di  complementaritd,  che  si  basa  sulfatto 
che  la  risoluzione  di  problemi  di  complementarila  si  pub  ricondurre  alia  risoluzione  di 
sistemi  di  equazioni  non  lineari  mediante  una  trasformazione  suggerita  da  Mangasar- 
ion.  R  metodo  e  “continuo’’  in  quanto  la  soluzione  del  sistema  viene  cercata  seguendo 
le  traiettorie  ottenute  per  integrazione  numerica  di  un'opportuna  equazione  differenziale 
—  come  in  precedents  lauori  degli  autori  —  e  Ji  pub  dire  “inesatto”  nel  senso  che  fa 
uso  di  un  metodo  di  gradienti  conitigati,  opportunamente  arvestato  “prima  della  con- 
vergenza”,  per  la  risoluzione  del  sistema  lineare  che  nasoe  nell 'integrazione  numerica 
dell 'equazione  differenziale.  R  metodo  appare  parlicolarmente  efficiente  per  problemi  in 
cui  compare  un  gran  numero  di  variabili  indipendenti,  net  quali  la  parte  prevalente  dello 
sforzo  di  calcolo  e  rappresentata  dalla  soluzione  di  un  sistema  lineare  ad  ogni  passo  di 
integrazione.  Vengono  dimostrate  la  convergenza  locale  e  la  convergenza  Q-superlineare 
del  metodo,  e  vengono  presentati  alcuni  risultaii  numerici  relativi  a  problemi  di  com¬ 
plementaritd  della  fisica  matematica. 

Abstract  -  We  consider  a  new  method  for  the  numerical  solution  both  of  non¬ 
linear  systems  of  equations  and  of  complementarity  problems,  based  on  the  fact  that 


^*^The  research  reported  in  this  document  has  been  made  possible  through  the  support 
and  sponsorship  of  the  U.S.  Government  through  its  European  Research  Office  of  the 
U.S.  Army  under  contract  n.  DAJA  4S-86-C-0028. 
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F.  ALUFFI-PENTINI  -  V.  PARISI  -  F.  ZIRJLLI 


[2] 


the  solution  of  complementarity  problems  can  be  reduced  to  the  solution  of  nonlinear 
systems  of  equations  by  means  of  a  transformation  first  suggested  by  Mangasarian.  The 
method  is  “continuous”  since  it  locks  for  a  solution  of  the  nonlinear  system  by  following 
the  numerical  solution  trajectories  of  a  suitable  differential  equation  —  as  in  previous 
work  by  the  present  authors  —  and  can  be  called  “inexact”  since  it  uses  a  conjugate- 
gradient  method  which  is  suitably  stopped  “before  convergence”  for  the  solution  of  the 
linear  system  arising  in  the  numerical  integration  of  the  differential  equation.  The 
method  appears  to  be  particularly  effective  for  problems  involving  a  large  number  of 
independent  variables,  where  the  computational  cost  is  dominated  by  the  solution  of  a 
linear  system  at  each  integration  step.  Local  convergence  and  Q-superlinear  convergence 
of  the  method  are  proved,  under  suitable  assumptions,  and  some  numerical  experience 
on  complementarity  problems  of  mathematical  physics  is  presented. 

Key  Words  -  Numerical  analysis  -  Nonlinear  equations  -  Mathematical  pro¬ 
gramming  -  Complementarity  problems. 

A.M.S.  Classification;  65H10  -  65K05 


1  -  Introduction 

Let  IR^  be  the  TV-dimensional  real  euclidean  space,  let  x  = 
(xi,X2,. ..  €  IR^  be  a  vector,  and  for  x,y  G  IR*'^  let  (x,y)  = 

N 

53  Xiyi,  |lx[|  =  (x,x)^'^  be  the  euclidean  scalar  product  and  norm;  where 

i=l 

necessary  ||  •  ||  will  indicate  also  the  matrix  norm  induced  by  the  eu¬ 
clidean  vector  norm.  Given  f  :  IR^  IR^  we  will  be  concerned  with  two 
classes  of  problems  in  this  paper:  the  problem  of  solving  the  system  of 
simultaneous  nonlinear  equations 

(1.1)  f(x)  =  0 

that  is:  find  x*  G  IR^  such  that  f(x*)  =  0,  and  the  complementarity 
problem 

(1.2)  X  >  0 

(1.3)  f(x)  >  0 

(1.4)  (x,f(x))  =  0 

where  x  >  0  means  Xj  >  0,  i  =  1,2, . . . ,  iV,  and  similarly  f(x)  >  0  means 
/i(x)  >  0,  »  =  1,2, . . . ,  iV,  /j(x)  being  the  components  of  f,  that  is:  find 
X*  such  that:  x"  >  0,  f(x*)  >  0,  (x*,f(x*))  =  0. 


[3] 


An  Inexact  Continuous  Method  for  the  Solution  etc. 
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The  importance  of  the  problem  of  solving  a  system  of  simultaneous 
equations  is  well  known.  When  f(x)  =  Ax  +  b  is  am  affine  map  the 
(linear)  complementarity  problem  has  been  considered  by  CoTTLE  and 
Dantzig  in  [1]  amd  contains  as  special  catses  the  lineatr  programming  and 
the  quadratic  programming  problem.  In  the  case  when  f(x)  is  a  possibly 
nonlinear  function  of  x  the  (nonlinear)  complementarity  problem  is  a 
rather  general  problem  and  contains  as  special  cases  the  Kuhn-Tucker 
first-order  necessary  conditions  for  the  nonlinear  programming  problem 
and  has  been  widely  studied;  see  for  example  Gould  and  Tolle  [2]. 

The  linear  and  nonlinear  complementarity  problems  have  applica¬ 
tions  in  such  diverse  areas  of  flow  in  porous  media  [3],  image  reconstruc¬ 
tion  [4],  [5],  game  theory  [6]. 

In  this  paper  we  will  be  concerned  with  the  problem  of  the  numerical 
solution  of  nonlinear  systems  of  equations  and  complementarity  problems. 
Usually  complementarity  problems  are  approached  numerically  with  piv¬ 
otal  methods  (for  example  the  simplex  method  for  linear  programming). 
The  pivotal  methods  are  usually  of  the  “step  by  step”  improvement  type, 
that  is,  given  a  problem  for  which  a  solution  is  sought,  the  standard 
approach  is  to  attempt  to  define  recursively  a  sequence  of  approximate 
solutions  which  have  the  basic  property  of  making  an  improvement  in  a 
suitable  “objective  function”.  When  the  problem  satisfies  some  convexity 
and/or  monotonicity  assumptions  the  pivotal  methods  are  guaranteed  to 
converge  and  if  only  a  moderate  number  of  independent  variable  is  in¬ 
volved  (up  to  few  hundreds)  their  numerical  performance  is  satisfactory. 

In  recent  years  there  has  been  a  growing  interest  in  the  use  of  con¬ 
tinuous  methods  in  nonlinear  optimization;  see  for  example  Allgower 
and  Georg  [7]  for  a  review  of  simplicial  methods  in  the  computation 
of  fixed  points  and  the  solution  of  nonlinear  equations,  and  Bayer  and 
Lagarias  [8]  for  the  interpretation  of  Karmarkar’s  linear  programming 
algorithm  as  a  method  that  follows  a  trajectory  of  a  suitable  system  of 
ordinary  differential  equations.  In  particular  the  present  authors  have 
developed  a  method  for  solving  systems  of  nonlinear  equations  based  on 
the  numerical  integration  of  an  initial-value  problem  for  a  system  of  or¬ 
dinary  differential  equations  inspired  by  classical  mechanics  [9],  [10],  [11], 
[12]  and  a  method  for  global  optimization  based  on  the  numerical  inte¬ 
gration  of  an  initial  value  problem  for  a  system  of  stochastic  differential 
equations  inspired  by  statistical  mechanics  [13],  [14],  [15].  In  section  2  the 
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algorithms  introduced  in  [10]  to  solve  systems  of  nonlinear  equations  are 
modified  to  obtain  an  “inexact”  solution  of  the  linear  systems  appearing 
in  each  iteration  in  the  spirit  of  Dembo,  Eisenstat  and  Steihaug  [16]. 
These  new  algorithms  are  particularly  effective  for  problems  involving  a 
large  number  of  independent  variables  where  the  computational  cost  is 
dominated  by  the  solution  of  the  linear  system  at  each  step.  Under  suit¬ 
able  hypotheses  local  convergence  and  Q-superlinear  convergence  of  these 
new  “inexact”  algorithm  for  nonlinear  systems  of  equations  dse  proved. 
In  section  3  the  complementarity  problem  is  transformed  into  a  nonlin¬ 
ear  system  of  equations  following  Mangasarian  [17]  and  therefore  the 
algorithms  previously  developed  provide  a  class  of  locally  convergent  Q- 
superlinear  methods,  which  are  not  of  the  “step-by-step  improvement” 
type,  for  the  solution  of  complementarity  problems.  Finally  in  section  4 
some  numerical  experience  obtained  with  the  algorithms  of  section  2  and 
3  on  some  complementarity  problems  of  mathematical  physics  is  shown. 

Some  of  the  results  of  this  paper  have  been  aimounced  in  [18]. 


2  -  Some  inexact  algorithms  for  nonlinear  systems  of  equations 

Let  f(x)=(/i(x),/3(x), . . .  ,//v(x))^e  IR^,  where  /,(x),  i=  1,2, . . . ,  iV, 
are  real-valued  regular  functions  defined  for  x  =  (zi,  zj, . . . ,  Xjv)^  6  IR^. 
In  order  to  solve  the  system  of  simultaneous  equations 

(2.1)  f(x)  =  0 
we  define 

(2.2)  f(x)  =  f(x)’-f(x)=f;/’(x). 

•  =  1 

It  is  easy  to  see  that  x*  is  an  isolated  minimizer  of  F(x)  and  F(x*)  = 

0. 

In  [9],  [10],  [11],  [12]  the  idea  has  been  proposed  and  developed  of 
associating  to  the  nonlinear  system  (2.1)  the  following  system  of  second- 
order  ordinary  differential  equations: 

(2.3)  M^(<)  =  -sD^W  -  Vf  (x(l))  I  e  (0,  +00) 
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Where  D  is  d.  N  x  N  positive  symmetric  matrix,  n,g  are  positive  con¬ 
stants,  VF{x)  is  the  gradient  of  the  function  F{x)  with  respect  to  x. 
The  equation  (2.3)  represents  Newton’s  second  law  (mass  x  acceleration 
=  force)  for  a  particle  of  mass  n  moving  in  subject  to  the  force  -Vf’ 
given  by  the  potential  F  and  to  the  dissipative  force  —gDdxfdi. 

If  x’  is  an  isolated  minimizer  of  F{x)  then  x{t)  =  x*,  V  t  E  [0,-|-oo), 
is  a  solution  of  (2.3);  consider  the  Cauchy  data 

(2.4)  ic(0)  =  fo 

(2.5)  §(0)  =  n, 

and  let  x(t,^o,t7o)  be  the  solution  of  the  Cauchy  problem  (2.3),  (2.4), 

(2.5) . 

It  can  be  shown  that  there  exists  a  neighborhood  U  C  IR.*‘'^  [o  ]  ^ 
such  that  if  E  f/  we  have: 

(2.6)  lim  llx(f,fo,  »7o)  -  x*l|  =  0 

*-•00 

Hence  in  order  to  solve  the  system  of  nonlinear  simultaneous  equa¬ 
tions  by  integrating  numerically  the  Cauchy  problem  (2.3),  (2.4),  (2.5), 
we  are  primarily  interested  in  the  equilibrium  points  reached  asymptoti¬ 
cally  by  the  trajectories  of  (2.3)  (hopefully  solutions  of  (2.1))  rather  than 
in  the  accuracy  of  the  numerical  scheme.  So  that  of  particular  interest  are 
numerical  methods  enjoying  a  special  stability  property  called  A-stability 
[10]. 

Let  t  E  IR,  let  y,fo  E  IR"*  and  <p{t,y)  E  IR"*  be  a  given  function 
continuous  in  t  and  continuously  differentiable  with  respect  to  y,  such 
that  the  initial-value  problem: 

(2.7)  ^(0  =  V>(«,y)  t€(0,-|-co) 

(2.8)  y(0)  =  eo 


has  a  solution  y(t,{o)  for  ^  €  [0,-1-co). 
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The  simplest  choice  of  A-stable  linearly  implicit  method  to  integrate 
numerically  (2.7),  (2.8)  is: 


(2.10) 


(/ - /i$„)(y„+i  - y„)  =  n  =  0,1,2,.., 

yo  =  fo 


where  y„  is  the  numerically  computed  approximation  of  y(nh,fo),  I  is 
the  identity  matrix  acting  on  IR”*,  h  >  0  is  the  stepsize,  for  n  =  0, 1,2, . . . 
u  =  n/i,  =  v{tn,yn),  =  Htn,yn)  wherc  $(t,y)  =  dipidy  is  the 
Jacobian  of  ip  with  respect  to  y.  We  note  that  when  9(t,y)  =  Ay  is  a 
linear  map  (2.9)  reduces  to  the  backward  Euler  method. 

After  rewriting  (2.3)  as  a  first-order  system 


(2.11) 


(2.12) 


^  -  -Vf’(x) 

dt  n  M 


formulae  (2.9),  (2.10)  with  variable  stepsize  h„,n  —  0,1,...  (i.e.  to  =  0, 
in  =  hi,  n  =  1,2,...)  are  applied  to  (2.11),  (2.12),  (2.4),  (2.5).  In  this 


case  the  map  <p:  IR^ 


IR^"  will  be  given  by 


(2.13) 


_ipv-iVF(x) 


so  that  its  Jacobian  matrix  is  given  by 


(2.14) 


$(x)  = 


0  I 

-H(x)  -iD 


where 


(2.15) 


I(x)  =  2  7(x)^y(x)  +  '^fi{x)ni{x) 


J{x)  =  df{x)/dx  is  the  Jacobian  of  f  with  respect  to  x  and  Hi{x)  is  the 
hessian  of  /i(x). 
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Let  s„  =  x„+i  —  Xn,  n  =  0,1,2,...;  after  some  simple  algebra  (2.9) 
becomes: 


(2.16) 


-Vf,  +  iv 
An 


n 


(2.17)  Vn+i  =  r-  n  =  0,1,2,... 

An 

(2.18)  x„+i  =  x„  +  s„ 

where  L„  =  L(Xr,},  Vf’„  =  VF(x„).  In  order  to  avoid  the  computation 

of  Hiix),  i  =  1,2, ...  ,iV,  at  each  iteration  and  since  we  are  looking  for 

s 

points  X*  such  that  f(x*)  =  0  the  term  /,(x)  in  (2.15)  is  dropped  so 
that  I(x)  is  substituted  by 

(2.20)  i(x)  =  2J’'(x)J(x). 

Equation  (2.16)  will  be  replaced  by 


(2.21) 

where 

(2.22) 

>l(x,/.)  =  i(x)  +  i[^/+!,D 

and 

(2.23) 

Afi  —  ^(XntAn) 

(2.24) 

I»„  =  ~VF„  +  -—Vn 

we  note  that  the  matrix  A„  is  symmetric  and  positive  definite. 
We  have  the  following  theorem: 
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Theorem  2.1.  Let  f  :  IR^ — be  twice  continuously  differen¬ 
tiable,  F{x)  =  f(x)^f(x)  and  L(x)  be  given  by  (2.15).  Let  x’  6  IR^  be 
such  that  f(x*)  =  0,  J(x*)  is  nonsingular  (i.e.  x‘  is  a  nondegenerate 
solution  of  the  system  (2.1)j  and  the  following  Lipschitz  conditions  holds: 

(2.25)  |ll(x)-l(x-)l|<l||x-x-|l  Vx6S  =  {x|||x-x-||<5} 

for  some  constants  7  and  6  greater  than  zero.  In  the  iteration  (2.21), 
(2.17),  (2.18)  let  {h„},n  =  0,1,2,...,  be  a  sequence  of  positive  numbers 
such  that 

(2.26)  lim  h„  —  00 

ft— *00 

then  there  exists  h  >  0  such  that  for  n  =  0, 1, . . . ,  x*  is  a  point 

of  attraction  of  {2.21),  (2.17),  (2.18)  and  the  rate  of  convergence  is 

(i)  Q-superlinear  if  h-'-  <  7itl^-f’(x„)|l,  7i  >  0,  n  >  Hq,  for  some  71, 
Tio  >  0- 

(ii)  Q-quadratic  if  h'^  <  72||Vf’(x„)jp,  72  >  0,  n  >  no,  for  some  72, 
no  >  0. 

Proof.  Let  us  rewrite  (2.21),  (2.17),  (2.18)  as 

(2.27)  x„^i  =  G(x„,  /in)  +  T—T  Af,  (x„  —  x„_i)  n  =  0, 1, 2, . . . 

nnnn-l 


where 

(2.28)  G(x,/i)  =  x-/i(x,/i)-'VF(x) 

with  the  initial  conditions  Xo  =  f 0  ~  ^-i  =  ^0, 

that  is  (2.21),  (2.17),  (2.18)  can  be  interpreted  as  a  two-step  iteration. 
Since  X*  is  a  nondegenerate  solution  of  the  system  (2.1)  x*  is  an  isolated 
minimizer  of  F{x)  and  VF{x')  =  0.  Moreover  for  /i  >  0  the  symmetric 
matrix  A{x,h)  is  positive  definite  so  that  A(x,/i)"'  e.xists  that  is  G(x,/i) 
is  well  defined  for  x  €  IR^  and  h  >  0  and  x’  is  fixed  point  of  G(x,  h). 
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Let  P  =  l(Z(x*)~^(l  and  let  e  €  (0,(2/3)"^)  then  there  exists  S  >  0 
and  h  >  0  such  that: 

|ll(x*)  -  A(x,/i)|l  <  £  Vx  €  S  =  (x|  llx  -  x'll  <  ^ j 

(2.29)  ^  ^ 

V/i  >  h 

In  fact 

lll(x-)  -  >l(x,/.)||  <  ||L(x-)  -  i(x)l|  +  ||i(x)-  A(x,ft)|| 
since  I(x*)  =  X(x*)  there  exists  S  such  that: 

||I(x-)  -  i{x)||  <  is  Vx£5 
and  for  a  suitable  h  >  0 

(2.30)  l|i(x)-4(x,/i)||=i||^/  +  sfll|<i;  V /.  >  S 

From  (2.29)  and  the  perturbation  lemma  (lemma  2.3.2  p.45  of  Or¬ 
tega  and  Rheinboldt  [19])  it  foUows  that  >l(x,h)"'  satisfies 

(2.31)  l|,i(x,/i)-‘|(  =  Vx65,V/»>h 

Moreover 

(2.32)  |lG(x,/i)-x*||  <u;(x,/i)|lx-x*l|  Vxe5,Vh>h 
where 

(2.33)  uix,h)  =  a[|M(x,A)  -  I(x)l|  -h  l|I(x)  -  Z(x-)|1  -h  ||9(x)!|] 
and 

fo 


9(x)  = 


||Vf(x)- VF(x-)-Z(x-)fx-x-)H 

lix-x'll 


X  =  X* 
X  7^  X* 
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In  fact 


llG(x,/.)  -  x'll  =  ||4(x,l.)-‘[4(x,/.)(x  -  X-)  -  Vf(x)]  11  < 

<  a{  [lM(x,ft)  -  I(x)ll  +  ||Z(x)  -  i(x-)||]  llx  -  x-||+ 
+  ||i(x-)(x  -  X-)  +  V/(x-)  -  Vf(x)||} 

Moreover  from  (2.25)  and  proposition  3.2.5  p.  70  of  [19]  we  have 


(2.34) 


lk(x)||  <  ai|lx-x*ll  Vx€5 


Hence  from  (2.30),  (2.25),  (2.33)  and  (2.34)  for  some  constants  03,03  >  0 
we  have 

(2.35)  ‘^(x,A)  <  oj^ oaljx  -  x’l)  Vx65,  V/»>h 

From  (2.27),  (2.31),  (2.32)  for  x,„Xn_i  €  S  and  >  h  we  have 

l|x„+i-x*||  < 

<  II G(x„ , )  -  X*  II  +  II  a; ‘[(x„  - X* )  +  (x*  - x„ _  1)]  II  < 

<  ["(*.,  k)  +  ^;^]  llx.. -X- II +^^lK..-x- 11  < 

<  a,«  +  a,i  +  ^j||x„-x'l|  +  ^||x,.,-x’|| 

Moreover  from  (2.36)  eventually  changing  the  values  of  6  and  h  we 


(2.37) 


so  that 


(2.38) 


.  02  uo  1 

7,  -  <..«  +  .J- + -p- <  5 

/iO  1 


Ix„+1  -  x*||  <  73||x„  -  x*||  +  74l!x„_i  -  x*| 


Ill] 
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with  =  73  +  74  <  1  that  is  x„+i  G  S.  In  particular  we  have  shown  that 


(2.39) 


lim  x„  =  X* 
n-^oo 


that  is  X*  is  a  point  of  attraction  of  (2.27). 

In  particular  for  n  >  no  >  0,  x„  6  5,  using  (2.36)  the  required 
order-of-convergence  estimates  follows  from: 


(2.40) 


||x„4.i-x’l|<  az^  +  oallxn -x’ll 


K  -x’||+ 


+ 


fia 


^n^n— 1 


|x„  -  x„_ijj  for  n  >  no  >  0 


and  the  fact  that 

(2.41)  l|Vf'(x„)ll<(i|i(x-)i|  +  s)llx„-x-ll 
where  lim  £„  =  0. 

n-*oo 

Using  the  method  given  by  (2.21),  (2.17),  (2.18)  requires  the  solution 
of  the  linear  system  (2.21)  at  each  step.  Computing  the  exact  solution 
with  a  direct  method  such  as  Gaussian  elimination  is  very  e.xpensive  when 
a  large  number  of  unknowns  is  involved  and  may  not  be  worthwhile  when 
X*  is  far  from  x* .  In  this  case  it  seems  natural  to  solve  the  linear  system 
(2.21)  by  an  iterative  procedure  and  to  accept  an  approximate  solution. 
In  particular  since  the  matrix  A„  is  symmetric  and  positive  definite  we 
may  use  conjugate  gradients.  When  the  method  given  by  (2.21),  (2.17), 
(2.18)  is  used  to  solve  (2.21)  with  an  iterative  procedure,  accepting  an  ap¬ 
proximate  solution,  we  will  describe  this  procedure  as  an  inexact  method. 

Let  s„  be  the  approximate  step  computed  by  the  iterative  procedure 
when  solving  (2.21)  and 

(2.42)  r„  =  AnSn  -  h„ 

be  the  residual.  When  r„  =  0  the  linear  system  is  solved  exactly.  Let 
us  assume  that  the  approximate  step  computed  s„  satisfies  the  following 
condition: 


(2.43) 


Iknll  <  PnWKW 


n  =  0,1,... 
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for  some  forcing  sequence  {Pn},  n  =  0, 1,. . .  We  have  the  following  theo¬ 
rem: 

Theorem  2.2.  Let  f  :  IR^ — be  twice  continuously  differ¬ 
entiable,  F{x.)  =  f(x)’'f(x)  and  £(x)  be  given  by  (2.15).  Let  x'  6  IR^ 
be  such  that  f(x*)  =  0,  J(x*)  is  nonsingular  and  the  following  Lipschitz 
condition  holds: 

(2.44)  lll(x)- I(x*)||  <  7|lx-x*ll  Vx  G  S  =  {x|llx-x*||  < 

for  some  constants  'f,6  greater  than  zero.  In  the  iteration  (2.21),  (2.17), 
(2.18)  let  {/in}i  n  =  0, 1,2,...,  6e  a  sequence  of  positive  numbers  and  let 
the  linear  system  (2.21)  be  solved  approximately  in  such  a  way  that  the 
residuals  r„  given  by  (2.42)  satisfy  the  condition  (2.43)  for  some  forcing 
sequence  {/?„},  n  =  0,1,....  If  0  n  =  0,1,...,  then 

there  exists  h  >  0  such  that  if  h„  >  h,  n  =  0,1,. . then  x’  is  a  point  of 
attraction  of  the  inexact  method  (2.21),  (2.17),  (2.18). 

Proof.  Since  /(x*)  is  nonsingular  and  L{x‘)  =  2J{x'ff  J{x')  we 
define  the  following  norm: 

(2.45)  11x11.  =  lli(x*)xll  VxGIR'' 
we  have 

(2.46)  ;7-IWl  <  ||x||.  <  ;„||x||  Vxem" 

Ml 

where 

(2.47)  p,=max{||i;(x-)||,||Z(x-)-'||} 

Moreover  it  is  easy  to  see  that  under  the  stated  hypotheses  for  any 
f  >  0  there  exists  ^  >  0  and  h  >  0  such  that: 

(2.48)  ll44(x,/i)- I(x-)11  <£  Vx€S=  {x|llx-x*|l<  tf}, /i>  A 

(2.49)  \\A{x,h)-^-L{xT^\\<€  VxG  S  =  {x|llx-x*ll  <6},h>h 
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(2.50) 


llVF(x)  -  VF(x*)  -  Lix)ix  -  x')||  <  £|lx  -  x*l| 

V  X  G  5  =  |x[  llx  -  x]!*  < 


We  have 
(2.51) 

i(x-)(x„+i  -  X-)  =  [/  +  L{x-){A-^  -  i(x*)-^)]  • 

•{r.  +  (>1„  -  L{x))iK  -  X-)  -  [  -  b„  -  VF(x*)  -  I(x*)(x„  -  x*)j } 

and  taking  norms; 


ll*n+l  -  X-||.  <  [l  +  ||L(X-)|1  IK*  -  i(xT 'll]  ■ 

(2-52)  ■[|lr4  +  p„-i(x-)ll||x„-x-|H- 

+11  -  b.  -  Vf(x-)  -  i(x-)(x„  -  x-)||| 
from  (2.24)  if  x„  G  S  and  h„  >  h  using  (2.48),  (2.49),  2.50)  we  have: 

llXn  +  l  -  x’l|.  <  [l  +  Mlc]  [/3r,llVF(x„)ll  +  2£l|x„  -  x‘||  + 


(2.53) 


I 


(1  +  /3n)(l|x„  -  X*||  +  ||X’  -  X„_i 


moreover  from 

(2.54)  VF(x„)  =  i(x*)(x„-x*)  +  [VF(x„)-VF(x*)-£(x-)(x„-x*)] 
we  have 


(2.55)  llVF(x„)ll  <  jlx„  -  x’ll.  +  £||x„  -  x*|| . 

Finally  from  (2.47),  (2.53),  (2.55)  we  have: 


|x„+i  -  x*||.[l  + 


/^max  (1  +  £/ii)  +  £^1(2  + 


/tVJ 


(2.56) 


lx„  —  X  ll.+[l  +  ;ii£]^^^(l  +  /?max)||x„_i  —  X* 


=  OollXn  -  x'll.  +  aeliXn-i  - 
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where 


(2.57) 


as  =  [1  +  Mi^] 


as  =  (1+/ii£)(1  +  /3 

m&x) 


choosing  the  values  of  e  and  h  so  that  as  +  as  <  1  from  (2.56)  we  have 
that  if  x„,x„_i  G  S  then  x„+i  €  S  and 


lim  x„  =  X* 

n— 


Theore.m  2.3.  Let  f  :  IR^ — be  twice  continuously  differ¬ 
entiable,  F{x)  =  f(x)”’^f(x)  and  L{x)  be  given  by  (2.15).  Let  x’  G 
be  such  that  f(x*)  =  0,  J{x')  is  nonsingular  and  the  following  Lipschitz 
condition  holds: 

(2.58)  ||i:(x)  -  i(x*)ll  <  7llx~x‘ll  Vx5  =  {x|llx-x*l|  <  ^} 

In  the  iteration  (2.21),  (2.17),  (2.18)  let  {h„],  n  =  0,1,...,  be  a 
sequence  of  positive  numbers  and  let  the  linear  system  (2.21)  be  solved 
approximately  in  such  a  way  that  residuals  r„  given  by  (2.42)  satisfy  the 
condition  (2.43)  for  some  forcing  sequence  {,3„,  n  =  0,1,...},  such  that 
0  <  <  /3tn«  <  1)  n  =  0,1, —  Then  there  exists  h  such  that  if  hn  >  h, 

n  =  0,1,..., X*  is  a  point  of  attraction  of  the  inexact  method  (2.21), 
(2.17),  (2.18)  and  the  rate  of  convergence  is: 

(i)  Q-superlinear  if  /i~‘  <  7i||VF(x„)||,  71  >  0,  n  >  no  for  some 
7i,  no  >  0  and  lim  ^„  =  0 

n— •« 

(ii)  Q-quadratic  if  h'^  <  72||V/’(x„)lp,  72  >  0,  n  >  no  and  /3„  < 
72llVF(x„)|l,  72  >  0,  n  >  no,  for  some  72,  no  >  0. 
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Proof.  From  Theorem  2.2  we  have  that  x*  is  a  point  of  attraction 
of  the  inexact  method  (2.21),  (2.17),  (2.18)  so  that  we  can  assume  that 
lim  x„  =  X*  and  it  remains  to  prove  the  rate-of-convergence  results. 

n-^oo 

We  have: 


(2.59) 


Xn+l  -X*  =  +  [a„  -  i(x*)](x„  -  x*)+ 

-  [  -  b„  -  VJ’(x*)  -  I(x*)(x„  -  X*)] } 


and  taking  norms 

ll*.«  -  x-ll  <  ||A;'||  [llr„||  +  II -  L(x-)||  ||x„  -  X-II+ 
(2-60)  +  ||Vf(x.)  -  Vf(x-)  -  i(x-)(x,  -  x-)||+ 


Let  e,6,h  be  chosen  in  such  a  way  that  (2.29),  (2.34),  (2.35)  hold 
then  there  exists  n'g  such  that  for  n  >  rig  +  1,  x„  6  5  =  {i]  ||x  -  x*||  < 
we  have: 


l|x„+i-x'll  <  a 


i3nllVF(x„)il+ 


(2.61)  +  ^02^  +  a3l|x„  -  x*l|^llx„  -  x'll  + 

+  Qi||x„  -  X*||^  +  ^(1  +  /3n)l|x„  -  X„_i|l 

and  the  desired  rate-of-convergence  results  follow  from  (2.41). 


3  —  Complementarity  problems  and  nonlinear  systems 

Let  f  :  IR^  — >  IR^  be  given,  the  complementarity  problem  cissoci- 
ated  with  f  is 


(3.1) 

(3.2) 

(3.3) 


X  >  0 
f(x)  >  0 
(x,f(x))  =  0 
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and  let  ^  :  IR  — *■  IR  be  a  strictly  increasing  function  such  that  ^(0)  =  0. 
In  [17]  Mangasarian  has  shown  that  x*  G  IR^  is  a  solution  of  the 
complementarity  problem  (3.1),  (3.2),  (3.3)  if  and  only  if  x‘  is  a  solution 
of  the  system  of  nonlinear  equations 

(3.4)  g(x)  =  0 
where  g(x)  =  (Si(x),52(x),...,5w(x))’'  and 

(3.5)  i?.(x)  =  <KI/.(x)-ar.|)-d(/;(x))-^x.)  i=  1,2,...,JV 
for  later  purposes  let  ns  introduce 

(3.6)  G(x)  =  g(x)’'g(x) 

DefIiNITION  3.1:  Let  x*  €  IR^  be  a  solution  of  the  complementarity 
problem  (3.1),  (3.2),  (3.3)  we  will  say  that  x*  is  nondegenerate  if  x’  + 
f(x*)  >  0. 

Definition  3.2:  Let  f  be  continuously  differentiable  and  J(x)  = 
df/dx  be  the  Jacobian  of  f  with  respect  to  x,  if  for  ii  =  1,2, . . . ,  iV  each 
principal  minor  {{dfi/dxj)),i,j  =  1,2,. ..,n,  is  nonsingular  we  say  that 
J{x)  has  nonsingular  principal  minors. 

In  [17]  Mangasarian  has  shown  that  if  x*  is  a  nondegenerate  solu¬ 
tion  of  the  complementarity  problem  (3.1),  (3.2),  (3.3)  such  that  J{x‘) 
has  nonsingular  principal  minors  and  0  :  IR  — ►  IR  is  a  strictly  increasing 
differentiable  function  such  that  d9ldt{0)  -|-  d9fdt{t)  >  0,  V  t  >  0,  then  x* 
is  a  solution  of  the  nonlinear  system  (3.4)  and  dg/dx{x’)  the  Jacobian  of 
g  with  respect  to  x  is  nonsingular. 

For  simplicity  we  choose  0(t)  =  t/2  so  that  in  a  neighborhood  of 
a  nondegenerate  solution  of  the  complementarity  problem  (3.1),  (3.2), 
(3.3)  the  function  g(x)  given  by  (3.5)  has  the  saime  regularity  properties 
of  f(x).  Given  the  local  character  of  the  convergence  theorems  of  section 
2  this  is  satisfactory.  In  section  4  the  method  for  solving  nonlinear  system 
described  in  section  2  will  be  applied  to  (3.4)  with  9{t)  =  t/2  for  some 
test  complementarity  problems. 
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4-  Numerical  experience 

The  inexact  method  (2.21),  (2.17),  (2.18)  has  been  implemented  as 
follows: 

i)  since  An  is  symmetric  and  positive  definite  the  linear  system  (2.21) 
has  been  solved  by  the  conjugate  gradient  method  (C.G.)  introduced 
by  Fletcher  and  Reeves  [20].  This  procedure  solves  a.n  N  x  N 
linear  system  in  at  most  N  steps.  Hovewer  we  stop  the  conjugate 
gradient  procedure  after  a  number  of  steps  which  is  usually  consider¬ 
ably  lower  than  N.  In  fact  let  be  the  approximate  value  for  the 
solution  s„  of  the  linear  system  (2.21)  obtained  as  the  result  of  step 
k  of  the  conjugate  gradient  iteration;  the  iteration  is  stopped  after 
step  m  if 

||.4„s<”>  -  b,||  <  ^„||b.|l 

ii)  We  have  chosen: 

fo  =  »?o  =  0 
p  =  P  =  1 

D  =  I  (the  identity’  matrix) 
sW  =  0  n  =  0,l,... 

and  the  following  very  simple  variation  laws  for  the  time  integration 
step-length  /i„  and  the  forcing  sequence 

hn  +  l  —  niin(10/ln, /iniax)  ^  ~  0,1,2,... 

with  /lo  =  1,  ~  10^* 

PUi=o„0l  n  =  0,1,2,... 

where  4o  is  given  and  d„  is  automatically  chosen  by  the  program 
among  the  two  .values  0.1  and  0.5. 

iii)  The  program  stops  in  any  case  the  conjugate-gradients  iteration  after 
N  steps  in  order  to  avoid  possible  non  termination  due  to  the  finite 
arithmetic  of  the  computer. 
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Finally  the  method  given  by  (2.21),  (2.17),  (2.18)  (i.e.  exact  solution 
of  the  linear  system  (2.21))  is  obtained  simply  setting  jSq  =  0- 

The  stopping  rule  adopted  is  C?(x„)  <  10“'°  for  the  inexact  method 
and  G(x„)  <  10“'°  for  the  “exact”  method  (i.e.  /3o  =  0).  These  methods 
have  been  coded  in  the  Pascal  programming  language  and  the  program 
has  been  run  on  a  Hewlett-Packard  9816  computer. 

We  have  tested  the  proposed  algorithm  on  three  complementarity 
problems  of  which  two  are  linear  and  one  is  nonlinear. 

The  first  problem  considered  arises  as  a  one-dimensional  free-boun¬ 
dary  problem  in  the  lubrication  theory  of  an  infinite  journal  bearing, 
i.e.  a  rotating  cylinder  separated  from  a  bezu-ing  surface  by  a  thin  film  of 
lubricating  fluid  [21].  The  finite-difference  approximation  used  by  Cryer 
in  [21]  leads  to 

Problem  A  (called  Problem  3D  by  Cryer);  Find  x,  w  g  IR^  such 

that 

(4.1)  w  =  q-f-A/x,  w>0,  x>0, 

(4.2)  (w,x)  =  0 

where  M  =  ((Af^)),  i,j  =  1,2, . . .  ,iV,  is  an  iV  x  A  matrix  with  elements 
Mij  given  by 


Afij  —  ~(^i+l/3)^  > 

if  j  =  i+l, 

if  j  =  i, 

Af,^  =  -(^,_i/j)", 

if  J  =  i  -  1 , 

O 

II 

otherwise 

and  q  =  (gi,  g:, . . .  ,qsY  is  a  vector  with  elements  g^  given  by 


(4.4) 

T  rr  ^ 

iV  1/21 1 

i  =  l,2,...jV 

where 

(4.5) 

T  1 

N  +  1. 
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and  the  function  n{y)  is  given  by 

(4.6}  ll{y)  = -^{1  + ecosTty)  >  0 

V'T 


with 

(4.7)  T  =  2,  £  =  0.8 

We  note  that  the  matrix  M  given  by  (4.3)  is  symmetric  and  positive- 
definite. 

The  second  problem  arises  as  a  two-dimensional  free-boundary  prob¬ 
lem  in  the  theory  of  the  steady-state  fluid  flow  through  porous  media. 
Some  of  these  problems  can  be  formulated  as  a  variational  inequality  af¬ 
ter  an  ingenious  transformation  proposed  by  Baiocchi  and  others  (ref. 
[13]).  The  discretization  used  on  the  “model  problem”  ([3],  p.  4)  leads  to 

Problem  B:  Find  x,  w  €  IR^  such  that 

(4.8)  w  =  q-|-Mx,  w>0,  x>0, 

(4.9)  (w,x)  =  0 

where  M,  an  iV  x  iV  real  matrix,  and  q  =  (gi,92»- ■  •  >9/^)^  €  IR^  are 
defined  below. 

Given  n,,n,  (positive  integers)  and  X,  Y  (positive  real  numbers),  let 

N  =  n^riy  , 

Dx  =  X/n,  -f  1 , 

Dy  =  Y/riy  -b  1 , 
a  =  Dy/Dx , 

let  A  be  the  n,  x  n*  tridiagonal  matrix  having  all  the  main  diagonal 
elements  equal  to  2(a  -f  1/a),  and  the  paradiagonal  elements  (i.e.  im¬ 
mediately  above  or  below  the  main  diagonal)  equal  to  —a,  and  let  B  be 
the  n,  X  n,  diagonal  matrix  with  diagonal  elements  equaJ  to  —1/a.  The 
matrix  M  is  an  n  x  n  matrix  with  a  block-tridiagonal  structure  (riy  x  n, 
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blocks),  having  each  main-diagonal  block  equal  to  the  matrix  A,  and  each 
paradiagonal  block  equal  to  the  matrix  B.  We  note  that  M  is  a  positive- 
definite  symmetric  matrix.  The  vector  q  is  defined  as  follows.  Given  W 
(0  <  W  <  Y),  and  using  the  Kronecker  symbol  6ij,  let 

gdy)  =  -  vY , 

9Riy)=\iW-y)\  if  y<W, 

gdy)  =  0 ,  if  y  >  w , 

5D(i)  =  yV2-(y^-w^)(x/2X), 

gu{x)  =  0 , 

Tij  =  -DxDy  -h  SiiagdjDy)  +  ^,n,ayn(;J?y)+ 

+  Sij{^l0')gD{iDx)  -f  6„^j{lla)gi;{iDx) , 

I  “  1?2,..., ,  j  —  1,2,..., Tty  . 

The  elements  gi ,  9: , . . . ,  of  q  a^e  given  by 

(4.10)  qk  =  Tij ,  with  k  =  {j  -  l)n,  +  i 

Our  last  problem,  which  is  defined  below,  can  be  interpreted  as  a 
finite-difference  approximation  of  a  nonlinear  variational  inequality. 

Problem  C:  Find  x,w  e  IR^  such  that 

(4.15)  w  =  Afx -t- p(x) -I- q ,  w>0,  x>0 

(4.16)  (w,x)  =  0 

The  problem  dimension  N,  the  quantities  Dx,  Dy  and  the  matrix 
M  are  defined  as  in  problem  B,  given  Ti,,n,,A’, y.  The  nonlinear  term 
p(x)  is  a  vector  in  IR^  with  components  pi  =  if,  :  =  The 

vector  q  =  (9i,9j, •  ■  •  ,9n)^  is  defined  by  equation  (4.10)  where  = 
DxDys\n{2TTiDxfX),  i  =  1,2, . . .  ,ti,,  j  =  1,2,.  ..,n,. 

The  numerical  results  obtained  with  the  previously  described  meth¬ 
ods  on  Problem  A,  B,  C  are  shown  in  Table  1.  2,  3  respectively. 
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Table  1  - 

Results  of  Problem  A 

no 

=  1 

no  = 

0 

n.  of 
step* 
(2.21) 

total  n. 
of  C.G. 
steps 

n.  of 
steps 
(2.21) 

total  n. 
of  C.G. 
steps 

30 

10 

79 

7 

210 

40 

12 

121 

8 

320 

SO 

16 

238 

8 

400 

60 

14 

240 

8 

480 

70 

IS 

318 

9 

630 

80 

IS 

369 

9 

720 

90 

19 

650 

9 

810 

100 

18 

SS6 

10 

1000 

Table  3  -  Results  of  Problem  B 
(with  X  =  1.62 ,  y  =  3.22 ,  ty  =  0.84) 


n. 

n. 

N 

no 

=  1 

no  =  0 

n.  of 
steps 
(2.21) 

total  n. 
of  C.G. 
steps 

n.  of 
steps 
(2.21) 

total  n. 
of  C.G. 
steps 

6 

9 

54 

13 

170 

6 

324 

8 

12 

96 

15 

250 

8 

768 

10 

15 

ISO 

17 

483 

10 

1500 

12 

18 

216 

19 

746 

12 

2592 

14 

21 

294 

19 

867 

14 

4116 

20 

30 

600 

34 

2405 

21 

126000 

Table  3  -  Results  of  Problem  C 
(withX  =  5.  y  =  5) 


n. 

"y 

N 

no 

n.  of 
steps 
(2.21) 

=  1 

total  n. 

of  C.G. 

steps 

no 

n.  of 
steps 
(2.21) 

=  0 

total  n. 
of  C.G. 
steps 

5 

5 

25 

5 

37 

4 

100 

10 

10 

100 

6 

99 

5 

500 

15 

IS 

225 

8 

278 

6 

1350 

20 

20 

400 

10 

407 

6 

2400 

2S 

35 

625 

10 

535 

8 

5000 

30 

30 

900 

10 

893 
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In  tables  1,  2,  3  the  adavantage  of  using  “inexact  linear  algebra”  with 
respect  to  complete  solution  of  the  linear  system  for  problems  A,  B,  C  is 
shown,  and  the  advantage  is  increasing  with  the  number  of  unknowns. 
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abstract 

A  new  method  to  solve  linear  programming  problems  is  introduced.  This  method 
follows  a  path  defined  by  a  system  of  o  d  e.,  and  for  nondegenerate  problems  is 
quadraticaliy  convergent. 


1.  INTRODLCTION 

Let  fl"  be  the  n-dimensional  real  Euclidean  space,  and 

x  =  (i| . x„V&R".  where  the  superscript  T  means  transpose.  For  x. 

y  e  fl"  let  x^y  be  the  usual  Euclidean  inner  product,  and  let  e  =(1 . 1)^. 
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A  linear  programming  problem  consists  in  minimizing  a  linear  function  oser 
a  region  defined  by  linear  equality  and  inequality  constraints.  We  will  say 
that  a  linear  programming  problem  is  in  canonical  form  when  it  is  written  as 
follows; 


minimize  c^x 
*e  «■ 

(1-1) 

subject  to 

Ax  =  0, 

(1.2) 

II 

(13) 

X  >  0, 

(14) 

with  side  conditions  Ae  =  0,  where 

.. 

rt  >  2 .  m  <  n. 

and  c  €  R"  are  given,  and  the  inequality  (1.4) 
that  is. 

is  understood  componentwise. 

It 

o' 

A\ 

.  11 . 

.Moreover,  we  assume  that  the  matri.x  .A  i>  of  rank  m. 

We  note  that  on  these  hypotheses  (l/n)e  is  a  feasible  point,  so  that  the 
feasible  region  is  not  empty. 

The  sintplex  method  applies  to  linear  programming  problems  in  standard 
form,  that  is. 

minimize  d’^v 

y  €  fi- 

(15) 

subject  to 

Cy  >  b. 

(16) 

y  >  0. 

(1.7) 

where  de  fi'.  befi',  Ce  fi'’*’.  r^s.  are  given,  and  the  inequalities  (1.6), 
(1.7)  are  understood  componentsvise.  In  [1]  it  has  been  shown  that  a  linear 
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programming  prohleni  in  standard  form  with  a  finite  solution  can  alwass  be 
reduced  to  canonical  form.  Moreover.  Kannarkar  assumes  that  the  objective 
function  of  the  problem  (1  l)-(1.4)  is  such  that 

a*  =  c’'x*  =  0 

for  any  feasible  point  x*  that  is  a  solution  of  the  linear  procramming  problem 
( 1.4).  The  problem  (1.1)-(1.4)  with  thi.v  extra  assumption  is  called  a 
problem  in  canonical  form  with  a  nonnalized  objective  function.  In  our  vxork 
TRe  assumption  of  havinc  a  normalized  objective  function  is  not  necessarv; 
however,  since  this  assumption  simplifies  some  of  the  following  alcebraic 
manipulations,  we  will  keep  it. 

Let  n  denote  the  siihspace  fl  ={xe  fl'M.Ax  =  0).  let  A  lie  the  simplex 
A  =  (x  e  fi"  ix  ^  0.  e^x  =  1).  and  rinaJly  let 

A  =  nnA  (IS) 

be  the  poiv  tope  of  the  feasible  points.  Then  the  problem  (l.l)-(l  4)  can  he 
rewritten  as  follows; 


minimize  c^x. 
1  e  \ 


(19) 


In  this  paper  we  will  introduce  a  new  method  to  solve  linear  program- 
ming  problems  in  canonical  form  with  a  normalized  objective  function.  This 
class  of  problems  is  the  one  considered  by  N.  Kamiarkar  in  his  celebrated 
paper  [2]. 

In  the  late  1940s  C.  B  Dantzig  [3]  developed  the  simplex  method  to 
solve  linear  programming  problems.  In  1972  V.  Klee  and  C  L.  Mints  [4] 
showed  that  the  worst  case  complexity  of  the  simplex  method  is  conibinato- 
rial.  Here  the  term  "complexity"  means  the  number  of  elementary  operations 
necessary  to  solve  a  linear  programming  pniblem  in  the  standard  form 
(I  5)-(1.7).  Since  the  simpic.x  method  finds  the  solution  after  finite  number 
of  iterations,  Klee  and  Minty  (4)  were  able  to  give  an  example  where  the 
simplex  method  has  complexity 


p  =  0{rs2'). 


Note  that  in  (1.5)  y  6  fi’.  Moreover,  in  1981  S.  Smale  in  (5]  showed  that  the 
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“average"  complexity  of  the  simplex  method  is 

p  =  0(rs*). 

where  r,s  are  the  dimensions  of  the  matrix  C  in  (1.6). 

In  spite  of  its  worst  case  combinatorial  complexity  ,  the  simplex  method 
has  been  very  successful  in  solving  linear  programming  problems.  The 
feature  of  the  simplex  method  responsible  for  its  worst  case  combinatorial 
complexity  is  that  it  moves  on  the  boundary  of  the  feasible  region 

Q  =  {y  e  fl'ICy  >  b,  y  >  0) . 

In  recent  >ears  a  great  deal  of  effort  has  been  spent  in  the  attempt  to  find  a 
new  algorithm  for  linear  programming  whose  complexity  in  the  worst  case  is 
polynomial.  It  is  believed  that  these  new  methods  will  go  through  the 
interior  of  the  feasible  region  Q. 

In  1979  L  C.  Khachijan  (6)  introduced  the  first  method  of  this  class, 
called  the  ellipsoid  metlwd.  The  worst  case  complexity  of  this  method  is 


P  =  0(s") 

Here,  however,  the  meaning  of  the  term  "complexity  '  has  been  slightly- 
changed.  In  fact  the  ellipsoid  method  does  not  step  after  a  finite  number  of 
iterations,  so  that  "complexity  means  the  number  of  elementary  operations 
necessary  to  arrive  in  a  predetermined  neighborhood  of  the  solution.  .More¬ 
over,  the  method  introduced  by  Khachijan  is  only  of  theoretical  interest, 
since  its  practical  perfonnance  is  rather  poor 

In  l98-i  .N.  Kannarkar  [2]  presented  a  new  linear  programming  method  of 
polynomial  worst  case  complexity 


p  =  0(s'’) 

This  algonthm  is  called  the  projective  method  when  applied  to  a  linear 
programming  problem  in  canonical  form  with  a  normalized  objective  func¬ 
tion.  This  algorithm  is  of  theoretical  and  practical  importance. 

Since  198-4  a  great  deal  of  work  has  been  done  in  developing  new- 
methods  for  linear  programming.  Several  "interior  point  algorithms'  have 
been  proposed.  P.  E.  Gill.  W.  Murray.  M.  .V  Saunders,  J.  A.  Tomlin,  and 
M.  H.  Wright  in  (7)  have  interpreted  Karmarkar  s  algorithm  as  a  "logarithmic 
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barrier  method"  and  have  suggested  a  new  algorithm  with  gcxxi  practical 
performance. 

In  the  framework  of  logarithmic  harrier  function  methods  we  can  recall 
the  work  of  several  authors.  In  [8]  J.  Renegar  lowered  Karmarkar's  comple.x- 
it>-  bound.  In  (9]  C.  Conzaga  lowered  Renegar's  comple.xity  bound.  In  [10]  M. 
Iri  and  H.  Imai.  with  the  hypothesis  of  being  able  to  perform  e.xact  line- 
searches,  introduced  a  Kiuadratically  convergent  algorithm  for  the  linear 
programming  problem.  In  [11]  N.  Megiddo  studied  the  geometr'''aI  proper¬ 
ties  of  the  paths  derived  from  "weighted  logarithmic  barrier  fut  'ticns." 
Finally.  J.  A.  Tomlin  in  [12]  reports  on  considerable  numerical  e.xperimenta- 
tion  with  this  kind  of  algorithms. 

In  this  paper,  as  suggested  by  D.  A.  Baser  and  J.  C.  Lacarias  in  [13].  we 
will  show  that  Karmarkar’s  projective  method  can  be  obtained  by  applying 
Eulers  method  with  variable  stepsize  to  a  suitable  initial  value  problen)  for  a 
system  of  ordinary  differential  e<|uations.  In  fact-  Karmarkar  s  method  obtains 
the  solution  x*  of  the  linear  programming  problem  b\  computing 

lim  x|f,  ( 110) 

where  x(t.(l/n)e)  is  the  solution  of  a  system  of  ordinary  differential  equa¬ 
tions  with  initial  condition  {l/n)e.  using  Euler’s  method  with  variable 
stepsize.  The  idea  of  obtaining  the  solution  of  nonlinear  programming 
problems  os  limit  points  of  the  trajectories  of  svstems  of  ordinary  differential 
equations  has  been  widely  used;  for  a  review  see  [14].  In  particular,  in 
[15-17]  quadratieally  convergent  algorithms  for  nonlinear  systems  of  equa¬ 
tions  have  been  obtained  from  methods  based  on  the  numerical  integration  of 
trajectories  of  systems  of  ordinary  differential  equations. 

The  interpretation  of  Karmarkar’s  projective  method  as  the  numerical 
solution  of  an  initial  value  problem  raises  two  natural  questions; 

(i)  Can  the  system  of  ordinary  differential  equations  used  in  Karmarkar  s 
projective  method  be  changed  to  a  new  one  that  will  generate  an  interesting 
algorithm’:* 

(ii)  Can  the  Euler  method  with  variable  stepsize  that  is  used  in  Kar¬ 
markar  s  projective  method  be  replaced  with  some  other  numerical  scheme 
that  will  generate  interesting  algorithms’ 

An  answer  to  question  (i)  has  been  given  by  D.  .A.  Bayer  and  J.  C. 
Lagarias  in  [13]  and  J.  L  Nazareth  in  [18],  who  replaced  Karmarkar’s  vector 
field  with  the  affine  teciof  field.  Question  (ii)  has  been  considered  by  N. 
Karmarkar,  J.  C.  Lagarias,  L.  Slutsman.  and  P.  W’ang  in  [19],  where  they  tried 
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to  approximate  the  p»ath  x(t.(l/n)e)  with  a  power  series  expansion,  obtaining; 
enctjuraijing  practical  results.  In  this  paper  we  give  two  new  answers  to 
(juestions  (i)  and  (iih  in  fact,  we  propose  a  vector  field  which  is  different 
from  the  ones  previously  considered,  and  we  use  a  linearly  implicit  .A-stable 
integration  scheme  [l-Jl  to  solve  the  initial  value  problem  considered.  In  this 
way  we  obtain  a  quadraticully  convergent  algorithm  for  linear  programming 
problem.  Moreover  our  algorithm  shows  good  practical  behavior. 

In  Section  2  Karmarkar’s  projective  method  is  interpreted  as  the  numeri¬ 
cal  integration  of  an  initial  value  problem  with  Euler’s  method  and  s-ariable 
stepsize.  Moreover,  to  a  linear  programming  problem  in  canonical  form  with 
normalized  objective  function  is  as.sociated  a  new  system  of  ordinarv  differ¬ 
ential  efjuations.  If  we  assume  that  the  solution  of  the  linear  programming 
problem  is  uni(]ue.  this  solution  can  be  obtained  as  the  limit  point  of  a 
suitable  trajectorx-  of  the  system  of  ordinary  differential  equations. 

In  Section  3  an  initial  value  problem  for  this  sy  stem  of  ordinary  differen¬ 
tial  etjuations  is  integrated  numerically,  using  a  linearly  implicit  .A-stahle 
method  with  variable  stepsize.  It  is  shown  that  this  is  a  (juadratically 
convergent  algorithm  for  linear  programming. 

Finally,  in  Section  4  we  compare  the  computational  cost  of  our  step  with 
that  of  Karmarkar  s  projective  algorithm  and  that  of  the  simplex  algorithm, 
and  we  present  some  numerical  experiments. 


2  THE  USE  OF  ORDINARY  DIFFERENTIAL  EQUATIONS 
(N  LINEAR  PROCR.3.M.M1.NC 

Let  X*  . t,;)'’  e  R".  fl'  Ive  given  by  A'*  =  (( .V,- ))  = 

Diagtx*).  that  is.  -A,*  =  i.J  =  1.2 . n.  where  S,j  is  the  Kroneckcr 

svmlvol. 

DtKisirios  2  1  .\  minimizer  .r*  of  the  problem  (1  1)-(1.4)  is  called 

nondegeiierate  if  it  has  exactly  »i  —  in  -  I  null  components. 

Let  y„={1.2. .  .  ii)  .ind  S  =  {s,.s. . s,.,  in 1  $  n.  be  an  ordered 

set  of  indices  such  that  S  C  Lc*t  z=(z,.;j . '  vector  U'e 

denote  bv  z.  the  vector  z.  =(c  . r  )6  fl'"*'.  Moreover,  given  a 

vector  v  S  ft'"*  ‘  and  a  matrix  of  rank  in  -  1.  we  denote  by 

the  submatrix  =  Iq  '.q’- . e  "  where  q'  is  the  jth 

column  of  Q.  If  B  is  an  ordered  set  of  indices  such  that  Bcj^  and 
.V  =  j„  —  B.  then  the  system 


Qz  =  V 


(2.1) 
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can  be  rewritten  in  the  following  fomj; 

Qb^b  +  Qs^s  =  '  (--) 

Defimtiov  2.2.  Let  S  be  an  ordered  set  of  m  ~  1  indices.  Then  B  is  a 
set  of  basic  indices  for  the  system  (2.2)  if  there  e.\ists  a  matri.v  € 
pdii-nxdi-m-i)  vector  vG  fi*"'  such  that  the  svstem  (2.2)  is  e()uiva- 
lent  to  the  system 

=  '  (-3) 

Lemm.-k  2.3.  Let  B  he  an  ordered  set  of  m  -r  \  indices  such  that  B  is  a  set 
of  basic  indices  for  the  system  (2.2).  Then  Qg  is  an  invertible  matrix. 

Proof,  it  follows  immediately  from  the  eciiiivalenee  of  the  linear  sv  stems 

(2.2)  and  (2.3).  ■ 

Ut  a>0.  x>0  x‘’=(.rf.x? . x^„f  e  B’'.  and  .V"  €  fi .  be  the 

matrix  A®  =  Diag(x“). 

Lemma  2.-1.  Let  a  >  1  arui  x*  be  a  nondegeneratc  minimizer  of  the 
linear  programming  problem  (!.!)-( 1.4).  Then  .AA’*“.A’'  is  an  invertible  ma¬ 
trix. 

Proof.  Let 


be  the  matri.x  .A  with  the  extra  row  added  Since  x*  is  a  nondegenerate 
minimizer  of  the  problem  (I.I)-(14)  and  M  has  rank  m  +  I,  then  there 
exists  an  ordered  set  of  m  1  indices  B  such  that  xj  e  fl'"*'  has  all  nonzero 
components  and  x*  €  fi""""'  is  the  zero  vector,  where  .V  =  ]„~  B.  More¬ 
over.  B  is  a  set  of  basic  indices  for  the  system 

Afy  =  u.  (2.4) 

where  uefi'"*'  is  given  and  y  G  R".  Ut  x-' e  fl'™ "  be  the 
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matrix  =  Diag<xj‘''-).  and  g  be  the  matrix 

X^'''*  =  Diag(x;.‘^^).  that  is,  the  null  matrix.  So  from  Lemma  2.3.  Mg  G 
R‘"* >>*("■*'>  is  invertible,  which  implies  that  SigXg'^^  is  invertible.  More¬ 
over,  since  B  is  a  set  of  basic  indices  for  the  system  (2.4),  we  have 


=  (.WgX;'/*  +  .v/v.v?'"")(x;'/-.vfj  +  x*'/*.v/0 

(2.5) 

Since  is  invertible,  from  (2.5)  it  follows  that  .V/XW/’’  is  invertible, 

so  that  an  easy  computation  shows  that  .AX*.A’’  is  invertible.  Therefore  it 
follows  that  .AX**''-  is  of  rank  m.  so  that  .A,V*“  ’  is  of  rank  m  and  .A.V*"A’' 
is  invertible.  * 


Lemm.a  2.5.  Let  a  =  l  or  a  =  2,  and  let  x*  be  a  nondegenerate 
minimizer  of  the  linear  programming  problem  (1.1)-(1.4).  Then  there  exists 
p*  >  0  sxich  that  .AX'.A’’  is  invertible  for  xGS(x*,p*).  where  S(x*.p*)  = 
{xefi'’lllx-x*||<p*}. 

Proof.  The  proof  follows  immediately  from  the  continuit\  of  .A.V^A'’ 
with  respect  to  x  €  fi'.  from  Lemma  2.4.  and  from  J.  M.  0rte2J  and  \\  C. 
Rheinboldt  [20.  Lemma  2.3.2.  p.  45).  * 

Let  X  =  (.t,.  t, . I,)'"  G  fi";  let  X  e  R"‘"  be  given  by  X  =  Diae(x);  let 

De  R'"*''.  m  ^  n.  bea  matrix,  and  D*  the  subspace 

D-  ={x  G  R-|Dx  =  0).  (2.6) 

and  let  )  be  the  orthogonal  projection  on  D".  The  projector  Ho-I  > 

alwavs  exists,  and  if  D  has  fiill  rank  is  given  by 

no.(.v)-[f-D''(DD^)’'D]y.  y  e  R".  (2  7) 


Let  r  be 


r  =  {xGRlx>0}. 


(2.3) 
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and  f  be  its  interior 

f-{iefi"l.x>0}.  (2.9) 

The  set  F  is  called  Ar  positive  orthant. 

For  a  >  0  and  ».€  F  let  A'"'*  e  fi"''"  he  the  matri.x  A’“  ’  = 
Diag(x“^-,jc“^-,....xJ-'*).  We  obsene  that  can  alwa\s  be  e.\- 

pressed  in  the  form  BT).  In  (act.  let  r  be  the  rank  of  the  matri.x  AX”  if 
r  =  m,  then  is  given  by  (2.7).  If  0  <  r  <  m.  we  can  consider  the 

matrix  AG  R'*"  obdiBed  from  A  by  eliminating  the  m  -  r  rows  of  .A  with 
indices  equal  to  thoir  of  the^m  — r  rows  of  .AV“'‘  that  are  linearly 
dependent.  Since  (.AX*  *)*  =(.A.V“'")*.  we  have 

“  ^^1  vv“  -)-■•■  (2.10) 

where  is  OMn  by  (2.7).  Finally  if  r  =  0  we  have  that  fluy.  =  /. 

where  /  is  the  n  X  «  idlentity  matrix.  Let  h(x)£  R"  he  the  following  vector 
field: 

h(x)*-A(/-ee^X)n,.,,v,-(Ac).  xSfl".  (2.11) 

The  vector  field  Wrf  is  known  as  Karmarkar  s  vector  field  [2.  13].  Let 
S(**.p*)  be  the  open  ball  of  Lemma  2.5,  and  let  us  consider  h(x)  for 
FuS(x*.p«). 

We  observe  that  br  x  €  FuS(x*.p*),  AA'  is  of  rank  m.  Then  h(x)  is  a 
continuously  djfferenlnble  function  of  x  for  x  €  F  U  S(x*,p*).  Let  A  be  given 
by  (1,8),  and  .\  be  tR«vn  by 

A  =  .\nf.  (2.12) 

W'e  will  consider  the  initial  value  problem 

dx 

■jj-  =  h(x).  (213) 

1 

x(0)=*-e.  (2.14) 

n 

It  is  easy  to  verify  tlnC  (l/n)e  e  A. 
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Lemm\  2.6.  Let 


B  = 


e  R‘*"* 


Dxn 


be  the  matrix  AX  uith  the  extra  rote  added,  and  let  y  =  x/^n{n-l). 
where  *  e  (0. 1)  is  a  parameter.  For  **  e  A  let 


=  Diag(x*). 


and  let  Xt^  be  given 


= 


(2.15) 


Then  Euler  s  method  applied  to  the  initial  value  problem  (2.13).  (2.14)  with 
variable  stepsize  Af*  given  by  (2.15)  produces  the  seiiuence  (x‘}.  k  = 

0.1.2 _  zenerated  by  Karmarkar's  algorithm  [2.  pp  3TS-3T9]  applied  to 

the  linear  programming  problem  with  normalized  objective  function 
(11)-(1.4) 


.  O 

Proof.  We  observe  that  At;  >  0  for  e  .\  (see  [2.  Theorem  5.  pp. 
351-382].  so  that,  integrating  (2  13).  (2.14)  with  Euler’s  method  and  variahle 
stepsize  At; .  we  have 


\ 


<) 


1 

—  e 
n 


x‘*' =  x‘ +  A(jh(x‘)  /t=0.1.2 . 


The  thesis  follows  from  a  straightforward  computation. 
For  a  >  0  let 


(2.16) 


(2  17) 


g(x.a)  =  n,^^.  :,-( .V"  '-c). 


xe  r 


(2.18) 
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Ue  note  that  for  a  =  1.  g<\.2)  ij.  the  affine  scaling  factor  of  [13].  and 

g(x.2)  =  n,,v,-(-Vc)  (2  19) 

can  be  defined  for  i  e  fi*  so  that 

h(*)  =  -  A•(/-ee^V)g(x.2).  xefl".  (2  20) 

We  have 

TiiEORtM  2  .  For  0  <  a  ^  2  let  X  ^  \  be  a  feasible  ptnnt  for  the  linear 
programming  problem  with  normalized  objective  function  (1.1)-(1.4).  Then 

n,.x- :,.(.V“  -c)  =0  (2  21) 

;/  arwf  oii’o  if  >  (j  n  iniiiintizer  of  the  linear  pro^ramminf’  problem  with 
normaitzed  objective  function  (1.1)-(1.4). 

Proof  Let  X  =  A'e  be  a  minimizer  of  the  problem  (1.1)-(14)  with 
normalized  objective  function.  Then 

n..o*^-(X“^-c)  =  0.  (2,22) 

In  fact,  if  we  assume  that 

n,,,- -c)  ^  0.  (2  23) 

then  there  exists  z  =  (z,.  . z,)€  fl"  such  that 

.  =  1 . m.  (2  24) 

;  -  I 

that  IS.  zeiAX'’  ')-  and 

m 

Zc,x:'-z,~0.  ^#0.  (2.25) 

/-I 


266 


S  HERZEL.  M  C.  RECCHIONI.  AND  F  ZIRILLI 


lhaf  is,  X“^^c  is  not  orthogonal  to  z.  We  can  assume  without 
generality  ^  >  0.  Let  us  define 

S  i  =  . n- 

Since  ^  >  0,  there  exists  j  such  that  *  0.  If  u.^  <  0  for  j  =  1.2. . 
choose  £  >  0;  otherwise  we  choose  c  as  follows; 

0  <  f  <  min  —  . 

j  u^>0  Wj 


We  recall  that  >  0  for  j  =  1.2 . n.  Let 

=  Xj  -  eWj.  j  =  1.2 . n. 


From  (2.2T)  we  have 


Vj>0.  j  =  1.2 . n. 


and 


I  S  >  0. 
;  - 1 


Let  us  define 


The  point  u  *(u,.u, . is  a  feasible  point;  in  fact. 

.Au  =  0. 

“  1. 

u  >  0 


loss  of 

(2.26) 

.,n.  we 

(2.27) 

(2.28) 

(2.29) 

(2.30) 

(2.31) 

(2  32) 

(233) 

(234) 
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Moreover. 


7-1 


e  V 


1 


e  V 


I  7-1 


=  — [  L  -£)3|  =  -}r-(c^X-  f^).  (2.35) 


e  V 


7  - 1 


e  V 


Since  x  has  been  assumed  to  be  a  minimizer  of  the  linear  programming 
problem  (l.l)-(1.4)  and  the  objective  function  is  normalized,  we  have 

c’’x  =  0.  (-2.36) 

Therefore  the  objective  function  assumes  a  negatise  value  at  u.  and  this  is 
absurd. 

Let  us  assume  now  that  x  G  .\  and  that  E(iualion  (2.21)  holds.  We  will 
show  that  X  =  .\'e  is  a  minimizer  of  the  linear  programming  problem  with 
normalized  objective  function  (l.l)-(l  4).  In  fact,  from  (2.21)  we  have 


(2.37) 

L'sing  A  instead  of  A  as  in  (2.10),  when  AX"'  ^  is  of  rank  less  than  m  we 
have 


0  =  e"X'-“'-n,,,.-  =,-(X“'-c) 

*  e’’X ' /  -  A"  'W( .A.\'"A’’)  ■ '  A.V“ X"  '=c 
=  e^Xc  -  e’‘XA’'(  .AX"A’’  )  ‘ '  .U'"c  *  e^’Xc  -  c’'x.  (2 .38) 


Therefore  we  have  that  x  is  a  feasible  point  where  c^x  =»  0.  that  is.  x  is  a 
•minimizer  for  the  linear  programming  problem  with  normalized  objective 
function  (l.l)-(l. 4).  ■ 
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Let  S  be  the  set  given  by 

S  =  {xefi«lAx  =  0.  e^-l}.  (2.39) 

Lemma  2.8.  Let  x*  be  a  nondegenerate  minimizer  of  the  linear  program¬ 
ming  problem  nith  normalized  objective  Junction  (1.1)-(1.4),  and  let  h(x)  be 
given  by  (2.11).  Then  we  have 

h(x*)=0.  (2.40) 

Moreover. 

Ah(x)  =  0.  xel.  (2.41) 

e'‘h(x)=0,  xSl.  (2  42) 

where  -  is  given  by  (2.39). 


Proof.  In  fact  for  x  €  S  we  have  .\Xe  =  0.  e^.Ve  =  1.  and 
.A.\nnv)-(.\'c)  =  0.  so  that 

.Ah(x)  =  -  .A.Vn,,^„.(  .Vc)  +  .-LVee^Vn,,,,.(  .Vc)  =  0  (2.43) 

and 

e''h(x)  =  -e^.\•^,.,„.(.\■c)+e^Xee^\•^,,,,.(■Vc)  =0  (2.44) 

.Moreover,  from  Theorem  2.7  we  have  (2  40).  ■ 

Let  X  e  fi".  and  £,  €  R"'"  be  the  matrix  given  bv 

£,  =  DiaE(n.,„.(.Vc)).  xefi".  (2  4.5) 

Let  /^(x)s  R’"”  be  the  following  matrix; 

y,(x)  =  -  ( /  -  XeeM.\  n,,,...\-'£,  -  [e'-.vn,,,,.(.xc)]/.  X  e  R\ 

(2.46) 
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For  X  e  R"  in  (2.46)  we  will  use  .A  instead  of  .A.  as  in  (2  10'.  when  .A.V  is  of 
rank  less  than  m.  An  elcinentars  computation  shows  that  the  niatri.\ 
.Vn,  .V  ‘  can  lie  defined  for  x  £  R"  so  that  /^(x)  is  defined  for  x  £  R". 
Lft 

■^  =  {x  £  /I'i-A-V'.A*^  is  invertible} .  (-47) 

For  x^  yf .  Ji,(x)  is  the  Jacobian  matrix  of  /i(x)  with  respctt  to  x.  .Moreover 
let  S(x*.p*)  be  tbe  open  ball  of  Lemma  2.5;  we  observe  that  for  x  €  F  U 
S(x*.p*).  since  .AA’  is  of  rank  in.  the  matri.v  .Vn,,y,-.V''  is  well  defined 

o' 

and  continuous.  So  J^Jix)  is  continuous  for  x  £  F  U  Six*,  p* ).  and  since 

FI,  =  0- 


7h(**)=0  (2.4S) 

From  Lemma  2.S.  we  conclude  that  any  solution  x*  of  the  linear  proviram- 
ming  prolilem  with  normalised  objective  function  (l.f)-(1.4)  is  an  eiiuilib* 
rium  pviint  of  (2. 1'}).  that  is.  h(x*)  =  0. 

However,  due  to  the  singular  Jacobian  of  h(x)  at  x*  [that  is.  to  (2.45)).  the 
use  of  a  linearlv  implicit  .A-stable  method  to  integrate  the  initial  value 
problem  (2.13).  (2.14).  as  suggested  in  (Sj  in  the  conte.vt  of  nonlinear 
programming,  will  not  produce  a  (juadratically  convergent  method  for  linear 
programming.  To  overcome  this  difficultv  we  introduce  a  new  vector  field 
Rxle  R"  defined  for  x  £  fl"  given  by 

f(x)= -(/-Xee'')[xc-XA’'(.A.\:AM’'.A.Vcj.  x  £  fi".  (2.49) 

"here  we  use  A  instead  of  .A.  as  in  (2.  10).  if  .A.X.a’’  .a  is  of  rank  less  than  m. 

Let  us  consider  fix)  for  x  £  F  U  S(x*. p*).  N\e  observe  that  .A.X.A^  is 
invertible  From  (2.49)  we  have  that  fix)  is  a  continuousK  differentiable 
function  of  x  for  x  e  F  U  S(x*.p*).  For  later  purposes  we  observe  that  fix)  for 
*  S  F  can  be  rewritten  as  follows; 

f(x)-  -(/-Xee'’)X'^-n,,v'  M-(-V'-c).  x  e  F.  (2.50) 
Or 

f(x)  = -X'  '-(/-X'-  'ee^X'''-)g(x.l).  x  £  F.  (2.51) 
Equations  (2  20).  (2.51)  and  Tbi'orem  2.7  it  follows  that  if  x*  is  a 
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minimizer  of  the  linear  programming  problem,  then  Rx*)  =  h(x*)^=  0,  that  is, 
X*  is  an  equilibrium  point  of  the  vector  fields  h(x),  Rx).  For  x  e  A  the  vector 
field  Rx)  can  be  obtained  as  the  steepest  descent  vector  associated  to  the 
function  c^x  with  respect  to  a  particular  metric.  In  [13]  D.  A.  Bayer  and  J,  C. 
Lagarias  have  introduced  the  idea  of  looking  at  Karmarkar’s  vec‘or  field  h(x) 
in  terms  of  steepest  descent  directions. 

Let  Xq  be  a  feasible  point  of  the  linear  programming  problem  (1.1)-(1.4), 
and  Fq  be  the  affine  subspace 

f„  =  x„-r{ve/l").Av=0.  e^v  =  0}.  (2.52) 

Lemma  2.9.  The  vector  field  Rx)  given  by  (2.51)  is  the  steepest  descent 
vector  associated  to  the  objective  function  b(\)  =  c^x  of  the  linear  program¬ 
ming  problem  (!.!)-( 1.4)  restricted  to  with  respect  to  the  Riemannian 

metric  C(x)=  .V"'  =  Dtag(x“').  defined  on  the  positive  orthant  P.  where 
is  given  by  (2.52). 

Proof.  We  consider  the  following  transformation  for  x  €  F^,^  T. 

X  =G'''-(x)y  =  .V'-v.  (2.53) 

5\’e  ha\e 

b{\{y))  =(.V'  'ej'^y  (2  ;  I) 

and  Ffl  assumes  the  following  form; 

Fu  =  x„t{u€R";.A.Y'"-u  =0.e’'.V'''-u  =0).  (2.55) 

The  gradient  vector  of  i;(x(y))with  respect  to  y  is 

db 

~  =  X''-c.  (2.56) 


The  gradient  vector  db/d\  projected  on 


is  given  by 
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where  we  use  A  instead  of  A.  as  in  (2.10).  if  AXA^  is  of  rank  less  than  m. 
Since  .AA’e  =  0  from  (2.57).  using  (2.7)  we  have 

5  =  .V'^c-  A•'''-ee^Vc.  (2.5S) 

Since  e^X.A^{AXA^)" 'AAc  =  o.  we  have 

4  =  (/-X''’-ee^A''^-)n,,,.,  (2,59) 

FinalK.  appKing  (2.53)  to  %.  we  have  that  the  gradient  vector  is  given 
by 

5(  x)  =  .V '  -i  =  .V"-(  /  -  A'-ee^V '-)  H, .V '  'c) ,  (2.60) 

This  concludes  the  proof  ■ 

Let  .\  be  given  by  (1.3). 

LtviviA  2.10.  Let  X  e  .\.  and  r  be  the  rank  of  the  matrix  AXA^ .  Then 

*(X'^'c)  =  (/-X‘''ee"A-''-)n,,,.  'c).  (2.61) 

uthere  we  use  A  instead  of  A.  as  in  (2,10).  when  r  is  less  than  m. 
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Since  AXe  =  0  and  e^Xe  =  1.  we  have 


A/  = 

(.AXA")" 

0 

0^ 

1 

(2.63) 


.^n  elementar>'  computation  gives  us 


n 


A.V''- 


e'.V‘ 


-(X'^'-c)  =  X‘^-c-  X'^U''(.-LX:a’’)‘'.A.\'c-  X'^-ee^Vc.  (2.64) 


From  e‘XA‘  =  0  we  have 


q  =  X '  ''-ee'’.X.A''(  AXA’^  ) '  '.AVc  =  0 


(2.65) 


and 


n 


A  V' 


-(X'"-c)  =  n.^v-  M-(  -V  'c)  -  X'^-ee^V'-c-q.  (2.66) 


With  an  eas\  computation  from  (2.66)  we  obtain  (2.61).  Let  r  =  0  and 
0  €  fi"  be  the  null  matri.v.  We  ha\e  that 


O 


.e^V'- 
50  that  (2.61 )  holds. 


=  (e^.V  *)  and  (1,;-=/. 


LtMMs  2.11.  Li’t  X*  be  a'nunde^encrate  minimizer  of  the  linear  pro- 
gruinminp  iiroblcw  with  normalized  objective  function  (1.1)-(1.4);  let  fix)  be 
given  by  (2.49)  and  1  be  given  by  (2.39).  Then  ut  have 


f(x*)  =  0. 


.67 


moreover 


Af(x)  =  0.  xel. 

e'’f(x)  =  0.  xel. 


(2  65) 
(2.69) 
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Proof.  Let  a  =  1  From  Theorem  2.T  we  have 

n,...v'M-(V*''-c)  =  0.  (2.70) 

so  that 

f(x*)=0  (2.71) 

Let  X  e  S.  Since  .A,Ve  =  0  and  e^.Ve  =  1.  we  have 

.Af(x)=  -(.A-.A.VeeM[.Vc-.>:A’’(A.\:A’’)''.A.Vc]  =0  (2.72) 

and 

e''f(x)  =  -(e'’-e\\'ee'’)[.Yc-  .\:a’’(  .A.MAM  ' '.A.Vc|  =0.  (2.7.3) 

This  concludes  the  pnwf.  ■ 

Lemn(v  2.12.  Let  Rx)  be  ptcn  by  (2.49).  and  x,,  G  fl’  be  such  that 

.Ax„  =  0,  (2.74) 

e'’x„=l.  (2.75) 

Then  the  solution  xit)  of  the  initial  lalue  problem 

dx 

—  =  f(x),  (2.76) 

x(0)  =  *„  (2.77) 

Mtisfiof  the  constraints 

Ax(t)  =  0.  (2.78) 

e''x(f)-l  (2.79) 

oU  values  of  t  where  x(<)  u  defined 


274 


S.  HERZEU  M  C.  RECCHIONI.  AND  F.  ZIRILLI 


Proof.  From  Lemma  2.11  we  have 

dx 

A— -Af(x)=0.  (2.80) 


dx 

e^—^e^f(x)  =  0,  (2.81) 

so  that  the  thesis  follows  immediately  from  the  assumption  (2.74).  (2.7.5)  on 
Xg  and  the  fundamental  theorem  of  calculus.  ■ 

For  xeR'let  £e/?"*"be  the  matrix  gi\en  by 

E  =  Diig[A^.\XA^)''AXc].  xefi".  (2.82) 

and  C  €  fl'’*"  be  the  matrix  C=  Diag(c).  Let  y(x)  e  fi"*"  be  the  following 
matrix: 

y(x)  =  -[/-.\lA'‘(.m'’)''A-Xee''](C-  £)^(e’’.X-c)/.  xefi". 

(2.8.3) 

where  we  use  .A  instead  of  .A.  as  in  (2.10).  if  A-\LA^  is  of  rank  less  than  m. 
Let  ^  he  the  set  (2.47).  .An  elementarv  computation  shows  that  for  xB 
y(x)  is  the  Jacobian  matri.x  of  Rx)  with  respect  to  x.  Moreover  let  S(x*.p*)  be 
the  open  ball  of  Lemma  2.5.  We  observe  that  for  x  e  F  U  S(x*.  p*).  since  the 
matrix  .A.X.A^  is  insertible,  J{x)  is  continuous  for  x  €  F  U  S(x*.p*). 

Theore.m  ;..13.  Let  us  assume  that  the  linear  procrammins  problem 
(1.1)-(1.4)  has  a  unique  nondegenerate  minimizer  x*.  and  let  J(x‘)  be  given 
,]y  (2.33).  Then  J{x*)  is  invertible  as  an  operator  restricted  to  the  subspace 

—  .  That  IS.  y(x*)v  *  0  for  each  v  #  0  such  that 


X 
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Proof.  First  of  al«r  show  that  for  i  =  1.2 _ .n  we  have 

( C— f)„  *  0  if  and  only  if  t  •  ^  0. 

Let  X*  —  Dia^x*).  Fob  Theorem  2.7  for  a  =  1  wc  have 

A-*'^-(C-f)-Diae((n.,v--  -•,-(A-''-c)))  =  0.  (2.S4) 

Since  C  -  E  is  a  diag^HJ  matrix,  from  (2.S4)  we  have  that  .Y  *  ^  0  implies 
(C— £),,  =0  for  i  *lL2....,n.  Let  us  show  that  (C-£),,  =  0  implies 

X*  *  0  for  i  =  1,2 - !».  In  fact  if  we  assume  that  there  exists  h  such  that 

(C  -  E)i,i,  =  0  and  then  from  the  assiinipticjn  that  x*  is  a  nondegen- 

erate  minimizer  of  thefcrar  programming  pmblem  (1.1)-(1.4)  it  follow  s  that 
there  exists  a  pivot  tasformation  that  makes  the  /ith  component  of  x* 
nonzero.  Let  y*  be  tl»-«rw  basic  feasible  solution  correspondinc  to  x*  via 
the  pivot  transformati»  Since  onl\  one  pivot  operation  has  been  made,  the 
nonbasic  components  Asr  than  the  hth  cx)mponent  are  still  nonbasic.  that 


is,  y*  =  0  for  each  i»lc  such  that  i*  =  0.  Let  V*  =  Diag(i/*,t/* . y*). 

From  (2,S4)  we  have 

y(C-£)=0  (2.&5) 


Moreover,  since  e^Vyf»0.  from  (2.85)  we  have 

0  =  e’^yfC  -  £)e  *i^r-Ce-e’’yA’'(  .A.y.A’' )  ' '.-LVc  =  e’’y*c  =  cNv 

(2.86) 

Therefore  y*  would  hranew  minimizer  of  the  linear  programming  problem 
(1  1)-{1.4).  different  x*.  and  this  is  absurd. 

We  have  /(x*)"  HC— £)".  In  fact  it  is  obvious  that  ve(C-  £)^ 
implies  v  e  J(x*)- .  Nlneaver.  let 

.\/-[4’'(.«'.A'’)''.A-ee^](C-  £).  (2  87) 

Then 

£)-X*.Vf,  (2  88) 

50  that  y(x*)v  =  0  im^^talC  —  E)v  =*  X*Mv.  Since  .Y*  is  a  diagonal  matrix. 
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we  obtain  (X*.V/v),  =  (X*)„(.Wv),,  i  =  1,2 . n.  We  have  two  cases; 

(i)  X*  =  0,  which  implies  ((C  -  £)v),  =  0; 

(ii)  X  •  *  0,  which  implies  (C  —  £)„  =  0. 


Summarizing,  we  have  that  {(C— £)v),  =  0,  i  =  1.2 . n.  that  is,  vS 

(C-  £)^. 

.Now  let  u  e  fi"  be  such  that 

.Au=0.  e’'u  =  0.  (2.S9) 


We  assume  that  u  6  J[x*)  ~ :  since  x*  e  (C  —  £)  ‘ ,  then  z  =  x*  +  u  G  ]{x*)  * . 
.Moreover. 

.Az  =  0.  e’’z  =  l.  (2  90) 

and  z  €  y(x*)  *  implies  z€(C  -  £)“.  If  ze(C  -  £)*.  then  =  0  for  each 
i  such  that  x  *  =  0;  this  condition,  together  with  (2.90).  is  a  characterization 
of  the  minimizer  x*  of  the  linear  programming  problem  (!.])-( 1.-1).  There¬ 
fore  u  =  0.  This  concludes  the  prcxif  ■ 

Theorem  2.14.  Let  x*  be  the  unique  nuru{e::cnerate  minimizer  of  the 
linear  programmmu  problem  uith  normalized  objective  function  (1.1)-{1.4). 
and  Rx)  be  given  by  (2.49).  ^^e  consider  the  initial  value  problem 

dx 

3-  =  f(x).  (2  91) 

at 


1 

x(0)  =  -e  (2  92) 

n 

Then  a  solution  .x(  t.  ( 1/ n  )e)  of  (2  91).  (2.92)  exists  for  t  €[0.  =).  and 


\  n  / 

Proof.  The  standard  existence  and  unirjueness  theorems  for  the  initial 
value  problem  for  ordinarv  differential  equations  guarantee  that  the  solution 
of  (2.91).  (2.92)  exists  locally  From  Lemma  2.9  it  follows  that  Rx)  is 
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tangential  to  S.\.  so  that  from  Lemma  2.12  and  the  fact  that  (l//t)e  E  A  ue 
have  that  x(f,(l/n)€)s  .\.  Moreover  for  x£  A  we  have 


d 

7t 


c^x 


-(.V  'c)  <0. 


(2.94) 


that  is.  the  objective  function  c^x  is  monotonicallv  decreasing  along  the 
trajectorv  x</.(l/n)e).  Since  the  minimunv  of  c^x  on  .\  is  zero,  x*  is  the 
unique  nunimizer  of  c^x  on  .\.  and  Rx*)  =  0.  from  C.  Sansone  and  R.  Conti 
[21,  p.  .3ll  vve  have  that  x*  is  the  unique  limit  point  of  x(/.(l/ /i)e)  and 


lim  X  t.  —  e j  =  X* 
(  -  *  »i  / 


(2.95) 


This  concludes  the  proof 
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J^t  xeR".  DcR"  be  an  open  set.  and  D  be  the  closure  of  D.  let 
^ .  D  C  R"  ■— R""  he  a  function  continuously  differentiable  in  D.  vvhose 
Jacobian  matrix  is  denoted  by  Q(x)  dw/dx.  Let  us  consider  the  initial 
value  problem 


dx 

■^  =  '^(x).  (3  1) 

x(0)  =  x„.  *„eD.  (3.2) 

Let  /  be  the  n  X  n  identitv  matrix,  >  0.  k  =  0.1.2 .  be  a  sequence  of 

stepsizes.  and  =  Then  any  solution  x(tj)  of  (3.1).  (3  2)  can  be 

Approximated  with  x*  computed  as  follows; 

x'’»x„.  (3  3) 

[/-/.i(?(x‘)ls‘»frtw(**).  t-0.1.2 .  (3.4) 


k  -0,1.2 . 


(3.5) 


27S 


S.  HERZEL.  M  C.  RECCHIONI.  AND  F.  ZIRILLI 


The  numencJ  scheme  (3.3)-(-3.5)  to  integrate  the  initial  value  problem  (3.1), 
(3.2)  is  .A-stable  and  linearly  implicit,  and  hasiseen  studied  by  J.  D.  Lambert 
and  S.  T.  Sigurdsson  in  [22]. 

Let  Rx)  be  the  vector  field  given  by  (2.49),  and  J(.x)  be  its  Jacobian 
matri-x  given  by  (2.83).  We  will  apply  the  numerical  scheme  (3.3)-(3.5)  to 
the  initial  value  problem 


dx 

3-  =  f(*). 
dt 

(3.6) 

1 

(0)  =  -e 

(3.T) 

n 


considered  in  Section  2.  Let  F  be  given  by  (2.9).  Bv  Lemma  2.5  there  exists 
a  neighborhaxl  S(x*.p*)  of  x*  such  that  fix)  is  a  continuously  diPTcrentiable 
function  in  F  o  S(x*.p*). 


Ltviviv  3  1.  Let  be  the  unique  tmnimtzer  of  the  linear  (programming 
problem  with  normalized  objective  function  (11)-(14).  Let  us  apply  the 
numerical  schernc  (3  3)-(3..5)  to  the  initial  value  problem  (3.6).  (3.7).  More¬ 
over  let 


hi^^{e^Xi^c)  for  0.1.2 . 

(3.S) 

1  ~  h^Jix'')  be  invertible  for  i=0.1.2 _ 

(3  9) 

Then  the  sequence  {x*  *'),  t  =  0. 1.2 .  generated  by  (3.3)-(3.5)  exists,  and 

$*  satisfies 

> 

l(  - 
o' 

II 

o 

(3  10) 

e’‘s‘=0.  it  =0  1,2 

(3  11) 

Proof  Ue  note  that  x"=(l/n)e  is  a  feasible  ptnnt 
programming  problem  (l.l)-(l  4)  and  that  s‘  is  defined  by 

of  the  linear 

{/-/.»y(x‘)]x‘-lij(x‘).  t=0.1.2 . 

(3.12) 

In  Lemma  2  11  it  has  been  shown  that  if  **  is  a  feasible  point  for  the  linear 
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programming  problem  (1.1)-(1.4).  we  have 


Af(x*)  =  0. 

-0. 

so  that  applying  A  to  both  sides  of  (3.12).  we  obtain 

A[/-hJ{*‘)]s*  =  0. 

From  (2. S3)  we  have 

.A7(*‘)  =  (e%c).A. 

Therefore  \^e  have 

(l-/it(e%c)lAs‘  =  0. 

so  that  from  (3.8)  we  have  As*  =  0.  1:  =  0, 1.2 . 

Moreover  it  is  easy  to  venfv  that 


e^;(*‘)=(e%c)e"/. 
so  that  from  Lemma  2.11  we  have 

e^[i-AJ(x‘)]s‘=0. 

"^^hich  implies 

(l-/.Je%c)]eV  =  0. 

so  that  from  (3.8)  we  have  e^s*  *  0.  it  *  0. 1.2 . 

Let  D  C  fl"  be  an  open  set  and  Dg  C  D  be  a  convex  set. 


(3.13) 

(3.14) 

(315) 

(3.16) 

(3.17) 

(3.18) 

(3.19) 

(3.20) 


Definitio.n  3.2.  Let  w(x);  D  £  R"  — •  fl*  be  a  continuously  differen- 
^ble  function.  Let  5  ^  R".  D^  be  an  open  neighborhood  of  4.  and  T .  Dx 
G  R"  X  R”  -»  L(R").  where  L(R")  is  the  set  of  the  n  x  n  matrices.  Then 
^^*■4)  is  a  consistent  approximation  to  the  Jacobian  matrix  Q(()  of  w(x)  on 
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D^C  D  if  0  e  fl"  is  a  limit  point  of  D^  and 

lim  r(x.|)  =  g(x)  (3.21) 

uniformly  for  x  e  Dg.  Moreover,  if  there  exist  t>*^t>  constants  c  >  0  and  r  >  0 
such  that 

||P(x)-r(x.4)||<c:li!  (3.22) 

for  each  x  e  D„  and  ^^D^nSlO.  r).  where  S(0.  r)  =  (^  e  fi"‘ |  ll|il  <  r). 
then  Tix.^)  is  a  strongly  consistent  approximation  to  Q{.\)  on  D„. 

Lem  MS  3.3.  Let  DQR"  be  an  open  set.  and  w:DcR''-*R"  be  a 
continuotisly  differentiable  function  on  the  conies  set  D,,  C  D.  Let  Q{x)S 
Lip^(  D„).  that  is.  let  Qix)  be  a  Upschitz  continuous  function  for  xE  with 
Lifischitz  constant  y  >  0.  Then 


]|  w( y)  -  w( x)  -  Q{ x)( y  -  x)  II  <  ^  be  -  y  il  ’ 


(3.23) 


for  each  x.  y  e  D,,. 


Proof. 
p  T3]. 


See  J  M.  Ortega  and  M  .  C.  Rheinboldt  [20.  Theorem  3  2  12, 


Theorem  3.4.  Let  x*  be  the  unique  nondezenerate  ininimizer  of  the 
linear  proitrainming:  problem  with  normalized  objectiie  function  (1  1)-(1  4). 
and  let 


B  = 


.A 


and  B  '  be  snen  by  t2  S).  Moreover,  let  fix)  be  ziien  by  (2.49),  and 


k  =0.1.2 . 


(3  24) 


where  Oj  >  o  >  0  is  a  bounded  sequence  such  that  h^  Then 
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there  exists  p,  >  0 open  neiahltorhood  S(x*. p, )  =  {*  e  fl "  i II*  -  x*ll  < 
Pi)  of  X*  such  th(it.Wli/£Six\Pf)<^  where  1  is  ^iien  by  (2.39).  then  the 
sequence  (x*),  h  generated  by 

(3.2.5) 

[/-/.iJtf§^-'-x‘)  =  /.J(xM.  Jt  =  0.1.2 .  (3.26) 

where  ]{\)  is  niicn^lBSSSi  am/  the  linear  system  (3.26)  is  suit  ed  in  B  ' .  is 

well  defined.  x^’’’^Cp,)ni  for  k  =  0.  1.2 .  and  {x^}.  k  =  0. 1.2 . 

contcraes  quadratiiM^  **. 

Proof  We  araJkinduefion  on  the  index  k.  Let 

I 

4-_.  L=u.l.2 .  (3.27) 

*4 

and 

<t>(i(^--fj^;(x‘)  il:=0.1.2 .  (3  2S) 

W’e  obserxe  that  (Skat  be  rewritten  as  follow  s: 

-<t>(^.0***’-x*)  =  f(xM.  A- =0.1. 2 .  (3.29) 

It  is  eas>  to  see  th^t®i^)  is  a  strongly  consistent  approximation  of  J(\)  on 
r  when  (  goes  toa*  We  haxe  seen  in  Lemma  3.1  that  if  v-e  fl"  and 
G  J  we  have 

JH**.f*)v={e"X.c-^i).Av.  (3.30) 

^»‘.fjv  =  {e%c-fje^.  (3.31) 

from  (3  24)  and  99iti]Baws  that  *  e^-X’^c  for  A  =  0. 1.2 .  so  that 

"'e  have 


if  and  onl\  if  v  e  B  *  . 


(3.32) 
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From  Theorem  2.13  we  have  that  J(x*)  restricted  to  B  is  invertible  and 

for  x*£S(z*,p,)  in  a  suitable  neighborhood  of  x*  and  in  a  suitable 
neighborhood  LI  of  zero.  The  perturbation  lemma  (see  J.  M.  Ortep  and 
W.  C.  Rheinboldt  (20,  Lemma  2.3  2,  p.  45])  implies  that  the  inverse  of  the 
linear  operator  <l>(x*,^j)  restricted  to  the  subspace  B*  e.xists  when  x*  e 
S(x*,p,),  €  U.  From  Lemma  2.11  we  have  that  x*  e  1  implies  Rx*)  e  B  *. 

When  X*  e  S(x*,p,)n  S  and  €  L’.  from  the  fact  that  Rx‘)eB'  and 
(3.32),  (3  33)  it  follows  that  x**‘  is  well  defined  and  x**'eS.  Moreover 
there  exists  tj  >  0  such  that 

|l‘t*(*-f)la*  |1< ’I-  xeS(x*.p,).  £eC’.  (3.34) 

and  we  have 

llx‘-'-x‘|l  =  ||4>(.x‘.^J|;.'[cb(x‘.^J(x‘-.t*)-f(x‘)||| 

<r,(l|D(x‘.f*)-/(x‘)|l  +  ||;(x‘)-/(x*);|)s‘-xM 

+  T,||f(.x*)-f(x‘)-;(x‘)(x‘-x*)||.  (3.35) 

Since  we  can  always  choose  Pj  >  0  such  that  y(x)e  Lip^(S(x*.p, ))  for  some 
y  >  0,  from  Lemma  3.3  we  have 


-x*ll  =Sf*77ll*‘  -**ll+  T7yllx‘  -x'l; 


+  ~'1*‘  -x" 


«a;(x‘,fj!!x‘ -x*!l.  (3  36) 

where 

inyllx‘  -x*ll.  (3.37) 

The  neighborhoods  S{x*.p,)  and  t’  can  be  chosen  in  such  a  wav  that 


"(**'f*)  <  *  <  1. 


x*6S{x*.p,)nS.  fjef.  (3.38) 
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From  (3.36)  and  (3. .36)  we  have 

-x*IUxllx‘-x«ll.  x‘ es{xvp,)ns.  (3,39) 

Therefore  x*'"'  e  S(x*. p, )r' I,,  Nioretner  we  ha\e 

!lx‘*'-x«IUx‘*'llx"-x*l|.  (3.40) 

so  that  the  iterates  {x*^},  h  =0.1.2 .  uiven  b\  (3.26)  are  well  defined  if 

x'’eS(x*.p|)nZ  and 

limx^=x*.  (3.41) 

I  -  X 

Moreoser.  since  Rx)  is  a  CDntiniiousU  di(TerentiaI)le  function  in  a  neichbor- 
hood  of  X*  and  fix')  =  0.  there  e.xists  a  constant  M  >  0  such  that 

||f(xM|i<  -.x*ll.  x‘ 6S(x*.p,)nl.  (3.42) 

From  (3.24)  we  obtain 

l|f(x‘)l|  M  ^ 

— ltx‘-x*l|.  x‘ eS(x*.p,)nI.  (3,43) 

w 

That  is.  the  sequence  (x*)  converges  quadratically  to  x*.  This  concludes  the 
proof  ■ 


T  NLMER1C.\L  E.XPERIMENTS 

N\e  begin  by  comparing  the  computational  cost  of  a  step  of  the  algorithm 
introduced  in  Section  3  with  that  of  a  step  of  the  simplex  algorithm  or  of  a 
step  of  Karmarkar’s  algorithm. 

consider  the  linear  programming  problem  in  the  canonical  form 
fl.l)-(1.4).  It  is  easy  to  verify-  that  the  computational  cost  of  a  step  of  the 
simplex  algorithm  is  given  by 


mn  +  lower  order  terms. 


(-1  1) 
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The  computatioittBSt  of  one  step  of  Karmarlur’s  algorithm  is  essentially 
due  to  the  computataaf  the  matrix 


AX\\^ 

and  to  the  solution iflle  m  linear  system 

(.■yf-A’’)y  =  .U'-c. 

Since  the  matrix  Ajar  is  symmetric,  its  computation  recjuires 

m~n 

— — y  lower  order  terms 


(-1.2) 


(4.3) 


(4.4) 


elementary  operat*  white  the  solution  of  the  linear  system  (4.3)  re<iuires 


m' 


— -  +  lower  order  terms 

6 


[4.5) 


elementary  operatin  Since  m  <  n.  we  can  conclude  that  the  computational 
cost  of  one  step  oflaBarkar's  algorithm  is  roughly 

f +  lower  order  terms  (4.6) 


elementary  operatic  This  cost  can  be  reduced  using  some  special  proce¬ 
dures;  for  exampleiB^  Karmarkar  has  shown  that  the  use  of  successive 
rank-one  modi ficat*«tt  compute  f4.2>,  (4.3)  reduces  the  "aserace  "  compu¬ 
tational  cost  of  eactep  to 


cn~^  *  lower  order  terms  (4.7) 

elementary  operatim  fcr  some  constant  c  >  0. 

The  computatiodtaBsC  of  one  step  of  the  algorithm  introduced  in  Section 
3  is  essentially  duetufce  computation  of  the  matrix 


.VUJ. 


(4.S) 
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to  thf  solution  ol  thi'  linear  s\  stem 

(.■L\:A'^)y  =  .-LVc.  (4.9) 

to  the  eomputation  of  the  niatri.v 

(4,10) 

that  appears  in  the  expression  of  the  Jatobian  Jix)  (2. S3),  and  to  the  solution 
of  the  linear  sxstem  {3.26).  The  computational  exists  of  (4.S).  (4.9)  are 
analo^uus  to  those  of  (4.2).  (4.3)  respectixely.  Moreover  the  computational 
cost  of  the  solution  of  the  linear  system  (3  26)  is 
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symmetric,  the  computational  cost  of  (4.10)  is 


m^n 


(4.16) 


elementary  operations.  Since  m  <  n,  we  cun  conclude  that  the  computational 
cost  of  one  step  of  the  algorithm  introduced  in  Section  3  is 

+  lower  order  terms  (4.17) 

elementary  operations.  Moreover,  if  we  use  successive  rank-one  modifica¬ 
tions,  as  proposed  in  [2],  we  can  decrease  the  "average  ”  computational  cost  of 
each  step  to 


c’n*  ’ -t  lower  order  terms  (4  IS) 

elementary  operations,  for  some  constant  c'>0.  Moreover,  to  improve  the 
value  of  c'  it  is  possible  to  use  anv  combination  of  the  ideas  proposed  in  [7. 
23.  24]. 

To  conclude,  the  computational  cost  of  one  step  of  the  algorithm  intro¬ 
duced  in  Section  >  is  of  the  same  order  as  that  of  one  step  of  Karmarkar  s 
algorithm,  vvhile  one  step  of  the  simplex  algorithm  is  much  cheaper.  How¬ 
ever,  due  to  ‘he  ■)l•.idratlc  convergence  of  our  algorithm,  we  expect  that  the 
numbt'"  of  jfe,  ..ns  needed  to  solve  a  linear  programming  problem  to  a 
given  accuracy  should  be  approximatelv  independent  of  the  problem  size  n. 

We  present  now  some  numerical  results  that  support  our  expectation 
The  algorithm  described  in  Section  3  has  been  implemented  using  two 
special  expedients  to  avoid  failure  due  to  the  ill-conditioning  of  the  problem 
considered. 

The  matrix  .AS  fl'"’"'  gj\(.n  by  (12)  is  replaced  with  the  matrix  .AS 
R"""  to  reduce  its  condition  number.  .A  is  obtained  using  the  singular  value 
decomposition  of  .A.A^  This  decomposition  has  a  computational  cost  of  order 
n’.  so  it  costs  the  same  os  one  step  of  the  algonthm  desenhed  in  Section  3 
Let  .A-A*^  =  Q^DV  be  the  singular  value  decomposition  of  .A.A^;  then  the 
matrix  A  is  given  bv 


(4  19) 

where  D*  G  fl"'*"'  jj  3  diagonal  matrix  such  that  (  D*),,  =  1/D,,  if  D,,  >  0 
and  (  D*)„  =  1  if  D„  =  0  for  i  =  1 . m. 
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In  the  first  k^  steps  of  our  algorithm  (i,  ^  5  in  our  numerical  experi¬ 
ments).  the  Riemannian  metnc 

C(x)  =  .V' =  Diag(*'')  (4.20) 

is  replaced  with 

G(x)  =  .\/,.V-'.  (4.21) 

where  .M,  =  .Xj''  =  Diagfx;"’)  and  X;  is  the  current  point  at  step  k.  We  note 
that  in  order  to  appl\  our  algorithm  is  not  nec-essar\'  to  have  a  normalized 
objective  function  —  that  is.  is  not  necessarv  to  know  the  value  of  the 
objective  function  at  the  minimizer  However,  in  our  numencal  e.vperin^ents 
we  use  test  problems  with  normalized  objective  function.  The  stopping  rule 
used  is 


t'i=c''x‘<1.0xl0''|c’’-J.  (4  22) 

W'e  have  considered  ten  test  problems  Problem  l(ZIR1)isa  problem 
coming  fpjm  the  operation  of  an  industnal  plant  in  central  Italy.  The  other 
problems  come  from  the  System  Optimization  Laboratorv  at  Stanford  Uni¬ 
versity  and  have  been  made  available  to  us  through  vtruiB  [25].  The 
numbers  of  variables  (n)  and  of  constraints  im).  shown  in  Table  1.  are  those 
relative  to  the  problems  in  canonical  form.  Finally,  it  denotes  the  index  of  the 
first  step  that  verifies  the  stopping  rule  (4-22). 

We  note  that  in  Table  1.  while  n.m  vary  by  an  order  of  magnitude,  the 
number  k  of  steps  needed  to  solve  the  problem  varies  only  by  a  factor  of  two. 
Moreover,  test  problems  with  n.m  ^  5  are  solved  in  about  ten  steps. 


TABLE  1 


Test  problem 

m 

n 

k 

Vi 

1. 

Z  IR1 

.304 

5-*3 

21 

2  41d-10 

2 

AOLI TTLE 

57 

141 

21 

3  16i>-09 

3 

AFIRO 

2S 

54 

12 

1  52t>12 

4. 

BEACONFO 

17.3 

298 

20 

3.91  d-09 

5 

BLEND 

75 

117 

21 

1.47d-12 

6. 

ISRAEL 

175 

319 

17 

1.94  D-IO 

1  . 

SC105 

106 

166 

13 

1.36i>ll 

S. 

SC50A 

51 

SI 

14 

1  48r>-14 

9 

SCS08 

51 

SI 

11 

7S4O-10 

10. 

SHARE2B 

97 

167 

21 

1  7bt>-10 
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The  algorithm  has  been  coded  in  fortran  and  tested  on  a  VAX/VMS 
Version  V5.1  in  double  precision  arithmetic. 
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